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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC) 
UNITS  OF  MEASUREMENT 


U.S.  customary  (non-SI)  units  of  measurement  used  in  this  report  can  be 
converted  to  metric  (SI)  units  as  follows: 


_ Multiply _ 

feet 

inches 

kip  (force)-ft 
kips  (force) 

kips  (force)  per  square  inch 
kip-inches 


_ 2z _ 

3.048 

2.54 

1355.818 

4.448222 

6.894757 


To  Obtain 
metres 
centimetres 
newton-metres 
kilonewtons 
megapascals 


112.9848 


newton-metres 


FOUNDATION  INTERACTION  PROBLEMS  INVOLVING  AN  ELASTIC  HALF-PLANE 


PART  I :  INTRODUCTION 

Background 

1.  Stress  analysis  problems  involving  interaction  between  a  structure 
and  its  foundation  often  lead  to  extensive  computational  effort.  This  fact 
becomes  obvious  when  finite-element  methods  are  used  to  study  problems  where 
the  foundation  is  idealized  as  an  infinite  elastic  half-plane.  Attempts  to 
represent  the  half-plane  by  a  finite  size  structure  using  many  elements  entail 
solving  large  systems  of  simultaneous  equations  at  a  considerable  computa¬ 
tional  cost. 

Purpose 

2.  This  report  presents  a  more  computationally  efficient  procedure 
which  employs  the  exact  solution  of  the  equations  of  elasticity  for  an  elastic 
half-plane  subjected  to  arbitrary  surface  loading.  Complex  variable  for¬ 
mulations  are  shown  to  yield  a  compact  solution  for  the  stresses  and  displace¬ 
ments  in  a  half-plane  supporting  several  concentrated  loads. 

3.  This  solution  is  also  employed  to  compute  flexibility  and  stiffness 
matrices  relating  the  concentrated  loads  and  the  displacements  at  the  points 
of  application  of  the  loads.  The  stiffness  matrix  derived  in  this  manner  is 
then  employed  to  investigate  a  simple  type  of  interaction  problem  where  an 

F.uler  beam  rests  on  an  elastic  half-plane  and  is  subjected  to  external  load- 

I 

ing.  The  foundation  interaction  forces,  as  well  as  other  quantities  such  as 
v'eam  shear  and  moment,  are  also  computed. 

i 


I 


PART  II:  TWO-DIMENSIONAL  ELASTOSTAT1C  PROBLEMS 


4.  A  brief  summary  of  the  various  field  quantities  important  in  the 
two-dimensional  infinitesimal  deformation  theory  of  linear  elasticity  is  pre¬ 
sented  below.  The  stress  state  is  represented  by  two  extensional  stresses 

and  t  and  a  shear  stress  t  .*  Hooke's  law  states  that  the 
xx  yv  xy 

stresses  depend  linearly  on  the  extensional  strains  £  and  t  and  the 

xx  yy 

shear  strain  £  .  Furthermore,  the  strains  are  functions  of  the  first 

xy 

derivatives  of  the  displacements  u  (in  the  x-direction)  and  v  (in  the 
y-direction) .  The  governing  differential  relations  are 


5.  Hook's  law  for  a  homogeneous,  isotrophic  material  contains  three 
elastic  constants  which  are  Young's  modulus,  E  ,  Poisson's  ratio,  o  ,  and 
the  shear  modulus,  G  .  Only  two  of  these  are  independent  because 


G  = 


E 

[2(1  +  o)] 


Since  a  two-dimensional  elasticity  problem  is  a  special  case  of  a  more  general 
three-dimensional  configuration,  it  is  instructive  to  examine  the  restriction 
involved  in  the  two-dimensional  specialization.  In  three  dimensions  we  have 


consider  the  following  quantities, 

each 

of  which 

depends 

on  spatial 

coordi 

tes  x  ,  v  ,  and  z  . 

a.  Three  displacements 

vi  , 

v  ,  w 

b.  Six  stresses 

T 

XX 

•  Tyy  ’ 

Tzz  ' 

T  ,  T 
xy  yz 

’  Tzx 

c.  Six  strains 

c 

XX 

ezz  * 

xy  yz 

'  %x 

+  S.  P.  Timoshenko  and  J.  N.  Goodier.  1970.  Theory  of  Elasticity,  3d  ed . , 
McGraw-Hill,  New  York. 


Subclasses — Plane  Strain  and  Stress 


fc.  The  main  subclasses  of  two-dimensional  problems,  plane  strain  and 
:e  stress,  are  mathematically  similar.  A  state  of  plane  strain  exists  when 
displacement  w  (in  the  z-direction)  vanishes  and  the  displacements  u 
v  are  functions  of  x  and  y  only.  It  is  then  found  that  for  plane 

tin 

e  =  c  =  e  =0 
yz  zx  zz 


T  =  I  =  0  ,  T  =  0(1  +  T  ) 

yz  zx  zz  xx  yv 


-cresses  T  ,  T  ,  and  T  related  to  strains  e  ,  E  ,  and  E 

xx  yy  xy  xx  yy  xy 

"ding  to 


2G  e  =  (1  -  o)x  -  T  ,  2G  e  =  (1  -  o)t  -  t 
xx  xx  yy  yy  yy  xx 


2G  f  =  x 
xy  xy 


.  tress  state  referred  to  as  plane  stress  occurs  when  all  field  quantities 

independent  of  z  and 


r  =  r  =  0  ,  t  =  0 
xz  yz  zz 


'  =  c  =  0  ,  E 

xz  yz  zz 


-0(T  +  T  ) 

_ xx  yy 


•  tinine  nonzero  field  quantities  T  ,  T  ,  T  ,  E 

xx  yy  xy  xx 

i  • e  related  for  the  case  of  plane  stress  by 


,  and 

yy 


2G  E 

xx 


-  OT 

X _ XI 

1  +  0 


2G  E 


yy 


T  OT 

yy  xx 

1  +  o 


2G  c  =  T 
xy  xy 
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7.  Because  the  equations  arising  in  plane  strain  and  plane  stress  are 
very  similar,  they  can  be  solved  by  the  same  mathematical  procedures.  It  is 
convenient  to  use  as  fundamental  elastic  constants  the  shear  modulus  G  and 
another  parameter  < 
where 


<  =  3  -  4o 


for  plane  strain  and 


< 


3-o 
1  +  o 


for  plane  stress.  Then  the  stress-strain  formulas  for  both  situations  turn 
out  to  be 


2G  c 


r  +  (<  -  3)  (t  +  t  ) 
xx _  xx  yy 


xx 


2G 


c  =  ^ 

yy 


+  (k  -  3)  (t 


XX 


2G  c  =  t 
xy  xy 

Although  the  constant  <  does  not  have  an  obvious  physical  significance  com- 
pirable  to  that  of  o  F  it  is,  nevertheless,  much  more  convenient  to  use  in  a 
varietv  of  formulas  needed  in  the  work  which  follows. 


Boundary  Value  Problems 


8.  Solving  a  boundary  value  problem  in  plane  elasticity  entails  deter- 


:  ng  the  stresses 


t  ,  and 
vv 


xv 


and  the  displacements  u  and  v 


result  from  given  boundary  conditions.  The  first  fundamental  type 


r<  b I  cm  involves  the  case  where  surface  traction  components  T  and  T  are 

xv 


These  tractions  depend  on  the  boundary  values  of  stress  according  to 


ve  n  . 


and 


T 

x 


T  V  +  T  V 

xx  x  xy  y 


T  =  x  v  +  i  v 
y  xy  x  yy  y 


with  and  being  the  x  and  y  components  of  the  outward  directed 

unit  surface  normal.  Assuming  that  the  boundary  tractions  comply  to  self- 
equilibrating  loads,  the  stresses  throughout  the  body  are  determined  uniquely 
and  the  displacements  are  determined  within  a  rigid  body  displacement  of  the 
form 


u  =  u0  “  V  *  v  =  vo  +  V 

where  u^  and  are  translation  components  and  a q  is  an  angle  of  rota¬ 

tion.  The  arbitrary  parameters  uq  »  vo  *  anc*  a0  can  ^xe<^  ^y  speci¬ 
fying  the  displacement  of  a  selected  point  and  the  direction  of  one  line 
segment  through  the  point. 

9.  Two  other  types  of  problems  also  merit  interest.  In  the  second 
fundamental  type  boundary  value  problem,  displacements  are  known  overall  on 
the  boundary.  Then,  the  stresses  and  displacements  will  be  uniquely  deter¬ 
mined  throughout  the  body.  A  third,  and  more  difficult,  type  problem  involves 
so-called  mixed  boundary  conditions  where  tractions  are  known  on  one  boundary 
part  and  displacement  conditions  are  known  on  another.  The  mixed  problem  is 
not  studied  here  extensively.  However,  one  important  mixed  problem  which  is 
dealt  with  to  some  extent  involves  a  half-plane  having  zero  load  on  one  part 
with  displacement  conditions  relating  to  an  Euler  beam  supported  by  the  half¬ 
plane.  To  analyze  this  problem,  developments  are  presented  in  the  following 
paragraphs  relating  to  a  half-plane  under  general  loading  and  an  Euler  beam 
under  general  loadings. 

Complex  Variable  Formulation 

10.  The  stresses  and  displacements  in  a  two-dimensional  linear  elasto- 
static  problem  can  be  expressed  in  terms  of  two  analytic  functions  <J>(z)  and 
f(z)  and  their  integrals*,** 


*  Timoshenko  and  Goodier,  op.  cit. 

**  N.  I.  Muskhel ishvili.  1953.  Some  Basic  Problems  of  the  Mathematical 
Theory  of  Elasticity,  Noordhoff,  Groningen. 


8 


$(z)  dz  and  iKz)  *  /t(z)  dz 


<J>  (z)  = 

When  the  boundary  conditions  and  the  boundary  geometry  have  simple  enough 
form,  general  solutions  for  $  and  ¥  can  be  written.  Included  in  this 
category  are  problems  for  the  circle  and  the  half-plane.  Several  stress 
formulas  needed  here  are  summarized  below.  These  and  other  such  relations  are 
developed  in  great  detail  in  the  reference  work  by  N.  I.  Muskhelishvili . * 

1 1 .  The  stresses  and  displacements  are  related  to  <f>  ,  f  ,  d>  ,  and 
^  ,  according  to  Kolosov's  formulas  which  are 


t  +  t  =  2 1  <J>(  z)  +  4>(z) 

XX  yy  L 

-t  +  t  +  2ir  =  2[zJ>'(z)  -+■  f(z)] 
xx  yy  xy 


2G(u  +  iv)  =  <<{> (z)  -  z$(z)  -  \p(z) 


where 


and 


i  =  /-i  ,  z  =  x  +  iy 


•c  =  3  -  4o 


for  plane  strain  or 


K 


3-o 
l  +  o 


for  plane  stress. 

12.  In  a  problem  leading  to  a  unique  stress  state,  then  4>  is  deter¬ 
mined  uniquely  except  for  an  additive  pure  imaginary  constant,  and  *1  is 


completely  determined.  Thus,  If  the  function  pairs  (♦,  'f)  and  ($q,  give 
the  same  stresses,  we  must  have 


4> 


0 


$  +  ic1 


V  =  t 
0 


where  is  real.  Integrating  for  and  f  shows  that 

♦  **  4>  +  ic^z  +  C£  ,  ^ 

with  c0  and  being  complex.  The  constants,  c^  ,  ,  and  c^  ,  corre¬ 

spond  to  a  rigid  body  displacement  field  of  the  form 


2G(u  +  iv)  =  (<  +  l)Cj (-y  +  ix)  +  (<c2  -  c^) 

In  a  problem  where  only  stresses  are  unique,  unique  displacements  can  be 
obtained  by  specifying  a  displacement  and  rotation  at  one  point.  Then  both  ♦ 
and  become  unique. 

13.  The  above  formulas  can  be  used  to  reduce  the  solution  of  plane 
elasticity  problems  to  an  equivalent  formulation  requiring  determination  of 
$(z)  and  T(z)  .  Furthermore,  these  functions  can  be  computed  in  general 
form  for  certain  restricted  types  of  geometries  and  boundary  conditions.  A 
case  of  special  interest  here  involves  the  elastic  half-plane  subjected  to 
stress-type  boundary  conditions.  Let  us  assume  that  a  half-plane  occupies  the 
region  -°°  <  x  <  °°  and  y  s  0  .  The  line  y  =  0  divides  the  plane  into  two 
parts : 

the  region  R+  defined  by  y  >  0  and 
the  region  R  defined  by  y  <  0 

The  stress  and  displacement  quantities  of  interest  will  exist  in  R  .  How¬ 
ever,  some  of  the  mathematical  formulas  presented  below  contain  quantities 
defined  for  both  R+  and  R  with  appropriate  limiting  values  being  used  for 
approaches  to  the  boundary  y  =  0  from  above  or  below. 

14.  Consider  the  case  with  the  normal  and  shear  stresses  known  at  all 


points  of  the  boundary.  Thus,  we  have 


.V.  aVa'a  \  OV  \’  - 


Tyy(t)  "  N(t)  Txy(t)  "  T(t) 


where  N  and  T  are  known  functions.  The  complex  stress  functions  which 
solve  this  problem  are 

00 

_  i  f  at  -  it)  dt 


1  f  (N  -  i1 

2iri  J  t  - 


00 

,/  x  1  f  (N  +  iT)  dt  A,  x  x 

'(z)  =  -  -tit  /  - — r - ' - ®(z)  ~  z«’(z) 

2ri  I  t  -  z 


Furthermore,  the  boundary  stresses  and  displacements  are  related  to  the 
boundary  values  of  $  in  a  remarkably  simple  form  according  to 

N  +  IT  «=  $+(t)  -  <f“(t) 

where  $  (t)  and  $  (t)  denote  limiting  values  of  ®(z)  at  a  boundary 
point  t  approached  from  above  or  below  the  x-axis,  respectively^ .  The  theo¬ 
retical  developments  relating  to  these  formulas  appear  in  Chapters  16  and  19 
of  Muskhelishvili's*  elasticity  book. 
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*  Muskhelishvili ,  op.  cit. 
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PART  III:  FUNDAMENTAL  LOAD  FORMULAS  FOR  THE  HALF-PLANE 


Stress-Function  Methodoloc 


15.  The  complex  stress-function  methodology  leads  to  simple  relations 
for  the  stresses  and  displacements  in  a  half-plane  subjected  to  certain  funda¬ 
mental  loading  conditions  such  as  application  of  a  load  distributed  over  a 
uniform  strip.  These  solutions  can  be  employed  to  develop  the  stiffness  and 
flexibility  matrices  relating  stresses  and  displacements  that  occur  when  a 
series  of  vertical  loads,  horizontal  loads,  and  couples  are  applied  to  the 
surface  of  a  half-plane. 

Vertical  and  Horizontal  Load  Solutions 

16.  Let  us  first  develop  a  solution  for  the  case  of  a  vertical  load  P 

and  a  horizontal  load  ,  distributed  uniformly  over  the  strip 


-a  <  x  <  a 


Thus,  we  have 


N  +  IT  =  0  for  I t I  >  a 


N  +  IT 


P0  +  iVo 


:  t  <  a 


v0  +  iP0 

$(z)  =  — -  [ &n (z  -  a)  -  Hn(z  +  a)] 

sra 


where  $(z)  is  an  analytical  function  in  the  plane  cut  along  a  straight  line 
between  z  =  -a  and  z  =  +a  .  The  formula  for  $  can  be  used  to  compute  the 
boundary  displacements  which  are  uniquely  determined  except  for  a  rigid  body 
displacement  and  rotation. 

17.  The  boundary  values  of  $  are  given  by 


where 


♦  "(t) 


a |  -  ln| t  +  a |  ±  i6(t) 


] 


Y0  = 


vo  +  iPo 

4tt  a 


and 


0(t)  =  0  for  i  1 1  >  a  and  9(t)  *=  *  for  |t|  <  a 


Consequently, 

2u  [u‘  (t)  +  iv'  (*>]  =  (<  +  1)Yq  [in  |t  -  a |  -  fcn|t  +  a|]  +  i(<  -  1)Yq" 


for  I t |  <  a 


and 


2u[u'  (t)  +  iv'  <‘>]  *  (K  +  l)Y^[lnlt  -  al-lnlt  +  alj  for  Itl  >  a 

18.  The  last  two  equations  can  be  combined  by  using  the  singularity 
function  <t  -  tg>n  which  equals  (t  -  tp)n  w^en  t  >  t^  and  is  zero 
otherwise.  Then  we  get 

u'(t)  +  iv'(t)  =  a^rCt)  +  i6^  |<t  +  a>0  -  <t  -  a>(^j 

r 

where 


(k  +  1)(V0  +  iPQ) 

°0  8u^a 


0  (K-  l>(Vo+iP0) 

60  "  2Pa 
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and 


r(t)  =  £n| t  -  a|  -  £n|t  +  a| 

Expressions  for  u(x)  +  iv(x)  can  be  obtained  by  integrating  u'  +  iv' 
is  not  hard  to  show  that 


J~  [<t  +  a>^  -  < t  -  a>^]  dt  =  <x  +  a>*  -  <x  -  a>*  -  a  =  g(x) 


/[■ 


£n|t  -  a  -  £n  t  +  a  dt  =  (x  -  a)  fcn  x  -  a|  -  (x  +  a)£n  x  +  a| 


+  2a£.n(a)  ■  f(x) 


so  the  general  displacements  equation  for  a  distributed  load  is 


u (x)  +  iv(x)  -  anf(x)  +  iBng(x) 


The  function  g(x)  is  odd  valued  and  has  three  different  intervals  of  defini¬ 
tion,  namely 


g (x)  =»  -a  , 


x  <  -a 


x  ,  -a  £  x  §  a 


x  >  a 


Furthermore,  f(x)  is  an  even-valued  function  having  a  logarithmic  singular¬ 
ity  at  | x |  *  ®  .  These  two  functions  form  the  basis  for  computing  the  defor¬ 
mation  effects  of  vertical  or  horizontal  forces  as  well  as  couples.  Several 
special  cases  arise. 


Case  I — Vertical  Force  Only  (positive  direction  upward).  In 
this  instance,  we  have 


V 


0 


0 


so  that 


iP0d  +  <) 

a  =  - - - 

0  8yira 


and 


$ 


0 


iP0(*  -  1) 
8pa 


Then 


u(x)  +  iv(x)  ■  i 


PQ(1  +  <) 
8u*a 


f(x)  - 


Vl  - » 

8ua 


g(x) 


The  real  and  imaginary  parts  of  this  equation  give 


Pn(l  +  O 

v(x) - sWT~  [f(x)  '  f(e)] 


and 


u  (x) 


~P0(K  ~  I} 
8Ua 


tg(*)l 


where  a  rigid  body  translation  term  has  been  added  into  the 
v(x)  equation  to  make  v  =  0  at  an  arbitrary  reference  point 
x  =  e  . 


Case  II — Horizontal  Force  Only  (positive  direction  to  the 
right) .  In  this  instance  P„  *  0  ,  so 


(k  +  1)V 


Ctn  = 


0 


(<  -  i)V 


0  8lJ7T 


and  = 


0 


0  8pa 


Then 


Vn(<  +  1) 

U(X)  -  -Tv—  [f(K)  -  f(6)] 


VQ(<  -  1) 

v(x)  =  -Sisi —  [8(x)] 


where  an  adjustment  has  been  made  in  the  horizontal  displace¬ 
ment  equation  to  make  u  vanish  at  x  =  e  . 

Case  III — Couple  (positive  direction  counterclockwise).  The 
load  effect  of  a  couple  can  be  obtained  by  combining  the 
effects  of  two  forces  of  equal  magnitude  and  opposite  direc¬ 
tion.  Let  us  represent  a  couple  of  magnitude  by  placing  a 

force  -M^/A  at  x  =  0  and  a  force  M^/A  at  x  *  A  .  We 

then  let  A  approach  zero.  The  resulting  displacement 
equation  is 


M0(<  +  U 

u(x)  =  — If(x  '  '  f(x)]/A 


which  becomes 


u(x) 


-Mn(iar)  -M-(ic  + 

“fcr  If,(x)]  'At  tr(x)1 


and  similarly. 


M0U 


1) 


Mo(<  r  0  0, 

-  [<x  +  a>  -  <x  -  a>  ) 


v(x)  = 


[g ’ (x)  ]  = 


Flexibility  and  Stiffness  Matrices 

19.  In  order  to  form  the  flexibility  and  stiffness  matrices  relating 
loads  and  displacements  for  a  series  of  vertical  forces,  horizontal  forces, 
and  couples,  it  is  necessary  to  know  the  deflection  equations  and  also  the 
rotation  equations  for  the  three  different  types  of  load  quantities.  It  is 
not  hard  to  see  that  the  rotation  of  any  initially  horizontal  element  on  the 
surface  of  the  half-plane  is 


e(x) 


dv  (x) 
dx 


so  the  rotation  function  is  obtained  by  differentiating  the  vertical  displace 
ment  equation. 

20.  The  purpose  of  distributing  forces  over  a  width  2a  is  to  avoid 
displacement  singularities  occurring  at  the  point  where  a  concentrated  load 
is  applied.  In  a  practical  application,  the  value  used  for  the  width  param¬ 
eter  should  be  small  compared  to  the  distance  between  successive  load  appli¬ 
cation  points.  Taking  a  =  d/20  ,  where  d  is  the  smallest  distance  between 
adjacent  load  points,  seems  to  be  a  reasonable  choice.  It  must  be  realized 
that  results  obtained  will  depend,  somewhat,  on  the  choice  of  parameter  a  . 
The  size  of  a  should  be  small  enough  so  that  a  concentrated  load  is  satis¬ 
factorily  replaced  by  a  distributed  load.  At  the  same  time,  this  parameter 
must  be  large  enough  so  that  round  off  errors  do  not  cause  inaccurate  calcula 
tion  of  the  various  influence  functions. 


Displacement  Loads 


21.  Let  us  now  summarize  the  various  displacement  formulas  correspond¬ 
ing  to  unit-load  quantities  applied  at  x  =  0  .  These  relations  involve 
elastic  constants 


c0 


<  +  1 

SUiS-  and  eo 


K  -  1 
8ua 


and  functions 


f(x)  =  (x  -  a)lnjx  -  a|  -  (x  +  a)£n|x  + 


r(x)  =  f'(x)  *  S.n|x  -  a|  -  fcn|x  + 


s(x)  =  f"(x) 


2a 


(x  -  a)  (x  +  a)  (x2  _  ^ 


g(x)  =  <x  +  a>*  -  <x  -  a>^  -  a 


h(x)  =  g'(x)  =  <x  +  a>^  -  <x  -  a>^ 


where  f  ,  s  ,  and  m  are  even  valued  and  r  and  g  are  odd  valued 
22.  For  a  vertical  load  at  x  =  0  we  have 


v  =  Cq [ f (x)  -  f (e) ]  =  CqF(x) 


u  =■  -eQg(x) 


9  --  cQ r(x) 


For  a  horizontal  load  at  x  =  0  we  have 


r’-f.  »• 


f.  ■' 


m 


mi 


ji 


.  <■  - 


v  =  eQg(x) 
u  =  CqF(x) 
9  =  eQh(x) 


For  a  couple  at  x  =  0  we  have 


I 

I 

i 

! 

i 


'A*  -  »v  •4,-'j.rv-.i,v.’  v.1' » 7-' 


23.  These  functions  can  be  used  to  formulate  the  flexibility  matrix  cor¬ 
responding  to  the  circumstance  where  horizontal  forces,  vertical  forces,  and 
couples  are  applied  at  a  series  of  points  on  the  surface  of  a  half-plane.  A 

program  is  given  below  which  forms  a  3n  by  3n  flexibility  matrix  for  a 

general  set  of  points  x^  .  x^  .  Horizontal,  vertical,  and  rotational 

quantities  in  the  matrix  correspond  to  row  and  column  positions  3i  -  2  , 

3 i  —  1  ,  and  3i  for  i  =  1  ,...,  n  . 


II*  ■  IMP  ■  Pf  V*L*\Vl*  W'.V 
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PART  IV:  ALTERNATIVE  REAL- VARIABLE  FORMULATION 
RESTRICTED  TO  VERTICAL  LOADS 


Complex  Variable  Methods 

24.  Complex  variable  methods  were  employed  above  to  derive  various  re¬ 
sults  applicable  for  a  half-plane.  When  only  vertical  loads  and  deflections 
are  of  interest,  then  equivalent  results  can  be  obtained  more  simply  by  using 
a  real-variable  formulation.  This  specialized  loading  condition  also  has  a 
meaningful  interpretation  in  a  three-dimensional  context  where  a  finite-width 
beam  rests  on  a  half-space  and  is  subjected  to  loads  symmetrical  about  the 
vertical  midplane  through  the  longitudinal  axis  of  the  beam.  The  vertical¬ 
loading  problem  will  now  be  studied  in  more  detail  before  moving  on  to  analy¬ 
sis  of  the  beam-interaction  problem. 

25.  Consider  an  infinite  half-space  the  surface  of  which  is  the  xz 

plane  and  the  interior  of  which  is  defined  by  y  <  0  .  We  assume  that  concen¬ 
trated  loads  in  the  y  direction  are  applied  at  positions  x^  , 

i  =  1  n  on  the  x-axis.  It  is  desired  to  compute  the  foundation  stiff¬ 

ness  matrix  K_  in  the  relation 
F 


v 


which  involves  the  applied  loads  P  ,  the  vertical  deflections  V  ,  and  the 

stiffness  matrix  K  .  Formulation  of  the  stiffness  matrix  is  to  be  based 
r 

upon  the  linear  theory  of  elasticity.  The  plane  analog  of  this  problem  occurs 
when  the  concentrated  point  loads  are  replaced  by  concentrated  line  loads  hav¬ 
ing  a  given  magnitude  per  unit  of  thickness  in  the  z-direction.  Each  of  these 
fundamental  problems  is  discussed  below. 

26.  The  deflection  functions  associated  with  a  concentrated  load  ap- 

*  **  f  ft 

plied  to  the  surface  of  a  half-space  are  well  known.  ’  ’  ’  When  a 


*  Timoshenko  and  Goodier,  op.  cit.  p  5. 

**  Muskhel ishvili,  op.  cit.  p  8. 

+  C.  S.  Desai  and  J.  T.  Christian.  1977.  Numerical  Methods  in  Geotechnical 
Engineering ,  McGraw-Hill,  New  York, 
ft  A.  P.  S.  Selvadurai.  1979.  Elastic  Analysis  of  Soil  Foundation 
Interaction ,  Elsevier,  New  York. 
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vertical  force  is  applied  at  some  surface  position  ,  then  the  verti¬ 

cal  deflection  occurring  at  another  surface  position  x_.  is  given  by 


a  -  a  )p0 
vij  “  "  V 


where  and  o  denote  Young's  modulus  and  Poisson's  ratio.  A  similar 

solution  exists  when  the  concentrated  point  load  is  replaced  by  a  concentrated 
line  load  distributed  parallel  to  the  z-axis.  In  that  instance  the  corre¬ 
sponding  formulas  relating  load  and  deflection  are 


2(1  -  a  )P, 


i.n|x^-  Xj j  for  plane  strain 


^  -  Xj 


for  plane  stress 


Infinite  Displacements 


27.  These  formulas  share  the  same  undesirable  characteristic  of  giving 
infinite  displacements  at  the  point  of  load  application.  This  difficulty  can 
be  remedied  by  averaging  the  load  over  a  finite  area  as  shown  in  Figure  1. 
Consider  the  case  where  a  concentrated  load  Pq  is  distributed  uniformly  over 
the  region  -a  ^  x  <  a  ,  -b<z<b.  Then  the  deflection  at  any  position  on 

the  x-axis  is  seen  by  superposition  to  be 


V1  - 

4rE„ab 


a 

ff 


dr  df 


2irE0ab 


J  {  y[77 

x-a  0  > 


Evaluation  of  the  above  double  integral  is  tedious  but  can  be  accomplished  in 
terms  of  elementary  functions.  It  turns  out  that 


b  +  (a 


(a  +  x)  .  b  +  \b  +  (a  +  x) 


b  a  +  x 


x  +  a  +\b^  +  (x  +  a)^ 


x-a  +\j b^  +  (a  -  x)^ 


In  the  above  formula  x  is  regarded  as  a  positive  quantity.  When  x  is 
negative,  |x|  should  be  used.  It  can  be  shown  that  as  x  becomes  large, 
the  last  formula  assumes  the  asymptotic  form 


v(x)  3 


V1  - 0  >  i 


as  would  be  expected  from  the  concentrated  load  solution.  In  fact,  when  x/a 
is  larger  than  six,  the  difference  between  the  concentrated  load  solution  and 
the  distributed  load  solution  is  negligible  regardless  of  the  value  of  b  . 

28.  A  distributed  load  solution  can  also  be  obtained  for  the  case  of 
plane  stress  or  plane  strain.  Let  us  assume  that  a  load  of  per  unit  of 

z-thickness  is  distributed  over  the  region  -a  <  x  <  a  .  Using  superposition 
for  the  plane  strain  solution  yields 


which  gives 


v(x) 


V1  -  °  > 

"Eoa 


x+a 
a 


V1  -  °  > 

7rE0a 


tf(x) 


t 

I 


2ai.n(a)  +  2a] 


where  f(x)  is  the  same  function  obtained  earlier  in  the  development  using 
complex  variables.  Hence,  the  displacement  formulas  given  by  the  two  methods 
are  the  same  with  the  exception  of  a  rigid  body  translation.  The  analogous 
formula  for  plane  stress  results  simply  by  replacing 


(1  -  o2)  with  1 


29.  The  distributed  load  solution  for  the  half-plane  is  bound  at  x  =  0 
but  has  a  logarithmic  singularity  at  x  =  «■>  .  As  a  remedy  for  this,  it  is 
reasonable  to  choose  a  value  e  which  is  much  larger  than  a  and  perform  a 
rigid  body  translation  such  that  v  =  0  at  x  =  ±e  .  Thus,  we  assume  that 
the  deflection  pattern  for  plane  strain  and  a  concentrated  line  load  is 


p0(1  -  0  > 

v(x)  =  - Tf— -  (f(x)  -  f  (e)  ] 

*Eoa 

where 

f (x)  =  -(x  +  a)S,n|  x  +  a|  +  (x  -  a)!,n|x  -  a| 

30.  In  the  instance  of  both  plane  loading  and  three-dimensional  load¬ 
ing,  the  formulas  for  deflection  at  a  distance  x  from  the  point  of  load  ap¬ 
plication  are  of  the  form 


v (x)  =  CqF (x  ,  a  ,  b) 

where  F  does  not  depend  on  the  elastic  constants.  The  constant  c  is 

o 

inversely  proportional  to  aE^  and  is  influenced  by  Poisson's  ratio  only  in 
the  cases  of  plane  strain  and  three-dimensional  loading.  By  superposition, 


the  surface  deflection  at  any  position  x  due  to  loads  P  at  x^  is 


v(x)  -  c0  ^  F (x  -  Xj  ,  a  ,  b)Pj 

j  =  l 


Evaluating  this  relation  at  x^  to  get  gives 


'i  ■  co  I  Vj 
3-1 


V  =  cqHP 


where  the  foundation  flexibility  matrix  H  has  elements 


hij  =  F^Xi  "  Xj  ’  a  *  b) 


Inverting  this  relation  yields  the  result 


P  «  K_V 
F 


where  the  desired  foundation  stiffness  matrix  K_  is  related  to  the  flexibil- 

F 

ity  coefficients  according  to 


K  =  —  H-1 
F  '0 


It  is  important  to  keep  in  mind  that  H  depends  only  on  dimension  parameters 

2 

a  ,  b  ,  and  ,  whereas  c^  involves  EQ  ,  (1  -  o  )  ,  and  parameter  a  , 

with  the  functional  dependence  which  these  variables  have  on  being  quite 

simple  in  form. 


PART  V:  ANALYSIS  OF  AN  EULER  BEAM  HAVING 
DISTRIBUTED  AND  CONCENTRATED  LOADS 


Solution  for  Load  Deflection  Relations 


31.  Before  an  analysis  can  be  made  of  the  interaction  between  a  beam 
and  a  half-plane,  necessary  load  deflection  relations  for  a  beam  are  needed. 
Toward  this  objective,  a  concise  solution  is  presented  below  for  the  shear, 
moment,  slope,  and  deflection  in  an  Euler  beam  having  constant  depth  and  sub¬ 
jected  to  an  arbitrary  combination  of  concentrated  loads  and  piecewise  lin¬ 
early  varying  (ramp)  loads.  This  solution  will  be  used  to  investigate  the 
interaction  between  a  loaded  beam  and  the  supporting  elastic  half-plane.  The 
contact  between  the  beam  and  the  half-plane  occurs  at  several  support  points 
which  can  transmit  only  vertical  concentrated  forces.  Displacement  continuity 
conditions  between  the  beam  and  the  plane  are  imposed  at  the  support  points. 

By  using  the  displacement  formulas  for  a  beam  and  a  half-plane  subjected  to 
the  same  concentrated  interaction  loads,  imposition  of  displacement  continuity 
conditions  at  the  interaction  points  gives  a  system  of  simultaneous  equations 
solvable  for  the  support  forces. 

32.  Beam  deflection  problems  can  be  conveniently  formulated  by  using 
the  singularity  function  <x  -  xQ>n  defined  such  that 

<x  -  xQ>n  =0  for  x  <  xQ 

=  (x  -  x0)n  for  X  g  Xq 

where  n  is  a  nonnegative  integer.  This  type  of  function  is  quite  useful  for 
describing  general  loading  conditions. 

33.  The  analysis  of  beam  deflection  problems  involves  the  quantities: 

load  per  unit  length  =  w 

shear  =  V 

moment  =  M 

slope  =  s 

vertical  deflection  =  v 

The  differential  equations  relating  these  variables  are 


where 


dV 

dx 


s 


y  L*'.*  W  *7*  !.■>  '.»*  I'T  V»  v»Hr 


dM  ds  M_  dv 

’  dx  “  V  '  dx  =  El  ’  dx 


E  =  Young's  modulus 
I  *  moment  of  inertia  of  the  beam 

Fundamental  Solutions 


34.  It  is  convenient  to  introduce  several  fundamental  solutions  used  in 
forming  more  general  solutions  by  superposition.  Consider  first  the  effects 
of  a  unit  load.  Assume  that  a  beam  occupies  the  region  x  ^  0  and  satisfies, 
at  x  =  0  ,  the  homogeneous  conditions  V=0,  M=0,  s=0,  and 
v  »  0  .  Also  assume  that  a  concentrated  load  of  unit  magnitude  acts  at 
x  =  c  .  It  is  found  that 


V  =  <x  -  c>^  ,  M  =  <x  -  c> * 

2  3 

l  <x  -  c>  =  l  <x  -  c> 

S  2  El  ’  V  =  6  El 

35.  The  load  and  deformation  quantities  associated  with  ramp  loading 
are  also  of  interest.  Assume  that  a  distributed  load  begins  at  x  =  a  with 
magnitude  R  and  varies  linearly  to  x  *  b  where  the  load  intensity  is  T  . 
The  slope  of  the  load  function  is  S  =  (T  -  R)/(b  -  a)  and  the  load  function 
for  this  case  is 

w  =  R  <x  -  a>^  +  S<x  -  a> *  -  T<x  -  b>^  -  S<x  -  b>* 


Integration  of  the  differential  equation,  subject  to  homogeneous  initial  con¬ 
ditions,  gives 


EIv  =  14 R<x  - a>A  +  no s<x  “ a>5  ~  k  T<x  ~ b>A '  no s<x " b>5 

36.  It  is  helpful  to  represent  the  previous  equations  in  the  following 
compact  form.  We  designate 

FUNCTION  FUNIT (X,C , El , ID) 

as  the  function  which  gives  the  effects  of  a  unit  concentrated  load.  The 
function  produces  shear,  moment,  slope,  or  deflection  according  to  whether 
ID  =  1  ,  2  ,  3  ,  or  4  ,  respectively.  Similarly,  we  designate 

FUNCTION  RAMP (X,A,R,B,T,EI,ID) 
as  the  function  returning  the  effects  of  a  ramp  loading. 

37.  The  fundamental  solutions  shown  above  can  be  used  to  formulate  and 
solve  a  general  problem  involving  multiple  loads.  Assume  that  the  following 
apply: 

a.  At  the  left  end  where  x  =  0  ,  the  moment  and  shear  vanish. 

b.  At  the  right  end  where  x  =  i.  ,  the  deflection  and  the  slope 
vanish . 

c.  At  several  points  external  concentrated  loads  act.  These  are 
described  by 

F^  at  dj  ,  j  •  1  ,  . . .  ,  nf 


At  several  points  foundation  reactions  occur. 


described  by  ,  ^  a(_  *  j  “  i  j  f  ii  •  me  iwuuuaiiuu 

reactions  are  initially  unknown  but  are  treated  as  parameters 
to  be  computed  later  using  displacement  constraints. 

Along  several  intervals  distributed  loads  act.  These  are 
described  by 


j  =  1  ,  ... 


These  are 
The  foundation 


Vx)  ■  Yx  -  */  +  Yx  -  ■/ 

-  T^<x  “  -  Sj<x  -  bj>*  »  j  =  1 . nr 


38.  Using  the  various  terms  developed  above,  the  deflection  equation 
for  a  beam  with  general  loading  and  general  end  conditions  is 


/  >  J  *  •  *  •  ' *  •  ’  •  ’  -  *  < *  "  .  •**  .  * 
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El  ,  4) 


vB(x)  =  vB(0)  +  Vg(0)x  +  vext(x)  +  2,  Pj *FUNIT (x  ,  Xj  , 

j-1 

where  Vg(0)  and  Vg(0)  represent  the  left  end  deflection  and  slope  which 
are  for  the  present  not  known  and 

nf 

VQV„(x)  =  Y  F  *FUNIT(x  ,  d,  ,  El  ,  4) 


PART  VI:  COMBINED  LOAD  DEFLECTION  RELATIONS 
FOR  THE  BEAM  AND  THE  FOUNDATION 


Analysis  Procedures 

39.  The  analysis  developed  in  Part  V  provides  the  tools  needed  to  form 
combined  flexibility  matrices  for  computing  interaction  forces  between  a 
loaded  beam  and  a  supporting  half-plane.  The  necessary  procedure  is  described 
in  the  following  paragraphs. 

40.  The  beam  is  subjected  to  a  general  external  loading  w(x)  plus  a 
series  of  reactions  P^  applied  at  the  foundation  support  points  x^  for 

i  =  1  ,  . . .  ,  n  .  The  boundary  conditions  on  the  beam  are  that  the  shear  and 
moment  vanish  at  each  end.  Integrating  the  differential  equation  for  the  beam 
deflection  gives  an  equation  of  the  form 


n 

vB(x)  =  vB(0)  +  v’(0)x  +  v£XT(x)  +  £  h.(x)P. 

J-l 

where  subscript  B  distinguishes  the  beam  and  P.  ,  ...  ,  P  ,  v„(0)  , 

1  n  B 

v'(0)  are  to  be  determined.  The  function  v_„_(x)  represents  the  combined 

D  tAl 

deflection  contribution  of  all  loads  except  for  the  foundation  reactions.  The 
function 


Y*> 


3 

<x  -  x .> 

_ J_ 

6EI 


represents  the  deflection  contribution  for  a  unit  load  applied  at  x.  .  The 
conditions  of  zero  shear  and  moment  at  the  right  end  are 


n 


I 


~VEXT(0 
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Where  VEXT 
loads . 


and  M, 


refer  to  shear  and  moment  terms  caused  by  the  external 


EXT 


Singularity  Avoided 

41.  Considering  the  foundation,  let  t(x)  represent  the  surface  de¬ 
flection  at  position  x  resulting  from  a  unit  load  applied  at  x  =  0  .  To 
avoid  the  singularity  in  the  concentrated  load  solution  for  a  half-plane,  it 
is  understood  that  the  concentrated  load  is  approximated  by  a  statically 
equivalent  uniform  pressure  applied  over  a  narrow  strip.  This  strip  is  con¬ 
sidered  narrow,  having  a  width  small  compared  to  the  total  length  of  the  beam. 
A  strip  width  of  11/1000  might  be  reasonable.  Using  superposition  indicates 

that  the  foundation  surface  deflection  caused  by  surface  loads  P,  ,  ...  ,  P 

I  n 

is 


j=n 

Vp(x)  =  -  t(x  -  Xj)Pj 

j=l 


Therefore,  the  flexibility  matrix  relating  surface  loads  and  deflections  and 
displacements  is 


w =  -  2 

j-1 


CijPj 


involving  a  symmetric  matrix  with  elements  t 


ij 


t(x1  -  Xj) 


A  negative 


sign  in  the  last  two  equations  is  used  because  an  upward  reaction  on  the  beam 
would  correspond  to  a  downward  (negative)  displacement  in  the  foundation. 
Matching  displacements  in  the  beam  and  half-plane  at  reaction  points  given 


t  •  •  •  9 


n 


^  t^ij  +  rij^Pj  +  VB^0^  +  VB^Xi  =  -VEXT^Xi^  ’  1  “  1  *  2 

j=! 


m 


where 


h 


ij 


h.(x.) 


is  supplemented  with  the  conditions  of  vanishing  shear  and  moment  at  the  right 
end,  a  system  of  n  +  2  equations  is  obtained  which  is  solvable  for  , 

P2  ’  Pn  *  VB^  ’  VB^  ‘ 

42.  Assuming  that  the  foundation  and  beam  load  transfer  mechanism  takes 
place  at  a  finite  number  of  support  points  amounts  to  an  effort  to  discretize 
a  continuous  foundation  pressure. 

Foundation  Pressure  Discretization 

43.  It  was  assumed  that  load  transfer  and  displacement  continuity  occur 
only  at  a  finite  number  of  support  points  along  the  beam.  In  this  manner,  a 
discretization  of  the  continuously  varying  foundation  pressure  is  achieved. 

For  satisfactory  results,  it  is  probably  necessary  to  use  40  or  more  interac¬ 
tion  points.  After  the  concentrated  load  reactions  have  been  found,  a  sta¬ 
tically  equivalent  piecewise  linear  load  distribution  can  be  derived  which 
better  represents  the  actual  interface  pressure.  With  all  loads  on  the  beam 
being  known,  any  additional  quantities,  such  as  internal  shear  and  moment,  can 
be  computed  for  arbitrary  positions  on  the  beam. 

44.  The  procedure  just  developed  employs  flexibility  matrices.  In  some 

cases  it  may  be  more  convenient  to  use  a  stiffness  matrix  formulation.  Such  a 

formulation  outline  follows.  Consider  a  beam  with  the  right  end  free  and 

left  end  values  of  shear,  moment,  slope,  and  deflection  which  are  F^  ,  -M^  , 

v'(0)  ,  and  v(0)  .  By  the  previously  discussed  methods  for  a  beam,  it  is 

found  that  the  deflection  v^  at  position  x^  associated  with  loads  , 

.  . .  ,  P  is 
n 
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=  v(0)  +  v'(0)xi  +  2^  bijPj 


where 


[xi  <3xj  -  xi> +  <xi  -  y3] 


Using  matrix  notation  we  can  write 


[vt  ,  ...  ,  vj1  =  V  ,  [v(0)  ,  v’(0)]T  =  VQ 


/  :  :::  ;  J  -  tF0  •  v  ■  p0 

I  n 


Assume  that  the  vector  of  nodal  forces  [P^  ,  ...  ,  P  ]  =  P  be  composed  of 

externa]  loads  P  and  foundation  reactions  P„  .  In  matrix  form,  the  beam 
fc-A  I  r 

deflection  equation  gives  V  =  RVQ  +  BP  which  also  implies  that 
P  =  K  V  -  K  RV  where  K  denotes  the  inverse  of  the  flexibility  matrix  B  . 

D  D  U  D 

Furthermore,  the  static  equilibrium  condition  balancing  the  left  end  loads  and 
the  external  loads  is  simply 

Po  .  -RTP  .  _8Tkbv  +  rTkbRV0 

Separating  the  load  vector  P  into  two  parts  gives 


?  =  K  V 

E  F  E  F 


where  denotes  the  foundation  stiffness  matrix.  Then  the  requirement  that 

the  foundation  and  the  beam  displacements  match  at  x,  ,  ...  ,  x  leads  to 
__  In 

Pr  -  K  V  =  K  V  -  K  RV  .  This  condition  coupled  with  the  requirement  that  P 


should  vanish  gives  the  following  system  of  equations  which  can  be  solved  for 

tv  ,  vy 


<kb  +  ke> 

-kbr" 

V 

V 

-<KbR)T 

T 

R \RJ 

vo 

0 

The  square  matrix  on  the  left  side  of  the  last  equation  is  symmetric  and  the 
system  can  be  solved  by  methods  such  as  Cholesky  decomposition  commonly  used 
in  finite-element  analysis.  Despite  the  fact  that  a  stiffness  matrix  approach 
can  be  used  in  analyzing  the  interaction  problem  of  interest  here,  using  flex¬ 
ibility  matrices  seems  more  natural  with  the  governing  differential  equations 
giving  deflections  directly  as  function  of  loads. 

Numerical  Results 


45.  The  analysis  procedures  developed  here  were  implemented  in  two  com¬ 
puter  programs.  The  first  program,  which  computes  shear,  moment,  deflection, 
and  interface  pressure  for  a  constant  depth  beam  supported  on  an  elastic  half¬ 
plane,  is  discussed  in  Appendix  A.  The  half-plane  can  be  in  plane  strain  or 
plane  stress.  The  beam  is  loaded  by  a  general  combination  of  concentrated 
loads  and  linearly  varying  distributed  loads.  The  interaction  between  the 
half-plane  and  the  beam  accounts  only  for  vertical  load  components  occurring 
at  an  arbitrary  number  of  supports  spaced  equally  between  limits  defined 
according  to  the  data  input.  A  system  of  equations  determining  the  interac¬ 
tion  forces  is  formulated  and  solved.  The  interaction  forces  are  then  re¬ 
placed  by  distributed  loads  to  provide  a  piecewise  constant  approximation  for 
the  foundation  pressure.  The  interface  pressure  is  used  to  compute  values  for 
the  shear,  moment,  and  deflection.  Three  example  problems  are  included  as 

typical  test  cases.  A  beam  1,200  in.*  long  and  120  in.  deep  has  a  modulus  of 

2 

5,000  kips/in.  and  rests  on  a  half-plane  having  an  elastic  modulus  equal  to 

2 

300  kips/in.  and  Poisson's  ratio  equal  to  0.2.  The  load  cases  studied 
inc lude : 


*  A  table  of  factors  for  converting 
ric)  units  is  presented  on  page  3. 


a.  A  500-kip  downward  load  is  distributed  uniformly  over  the  total 
length.  (The  beam  width  is  1  in.  measured  normal  to  the 
xy-plane . 

b.  Two  concentrated  downward  loads,  each  having  a  magnitude  of 
250  kips,  are  applied  at  the  beam  ends. 

£.  Two  concentrated  couples,  each  having  moment  magnitudes  of 

6,000  kip-in.,  are  applied  at  the  beam  ends.  The  sense  of  the 
couples  at  the  left  and  the  right  ends  are  clockwise  and  coun¬ 
terclockwise,  respectively.  The  couples  are  simulated  by  using 
ramp  loads. 

46.  Computer  output  for  these  problems  appears  in  Appendix  A.  These 
results  exhibit  the  type  of  behavior  typically  observed  in  punch  problems, 
namely  the  occurrence  of  very  large  stresses  near  the  ends  of  the  beam.  In 
the  case  of  the  uniformly  loaded  beam,  the  large  interface  pressure  at  the 
ends  rapidly  diminishes  to  a  nearly  uniform  value  near  the  middle  of  the  beam. 
When  the  same  external  load  resultant  is  concentrated  at  the  ends,  rather  than 
being  distributed  uniformly,  most  of  the  load  goes  into  the  foundation  at  the 
ends.  Because  the  half-plane  tends  to  deflect  downward  at  the  midpoint  more 
so  than  the  beam  does,  a  reversal  in  the  reaction  direction  occurs  near  the 
beam  center  to  maintain  displacement  continuity.  In  the  third  case,  the  two 
concentrated  end  couples  have  a  zero  statical  resultant  and  tend  to  cause  the 
beam  to  have  compression  in  the  top  and  tension  in  the  bottom.  If  the  beam 
ends  were  not  connected  to  the  half-plane,  the  ends  would  rise.  Consequently, 
imposition  of  displacement  continuity  causes  a  tensile  interface  reaction  at 
the  ends  and  compression  at  the  middle. 

47.  Some  discussion  of  the  approximations  used  to  obtain  distributed 
interface  pressures  should  be  presented.  After  the  reactions  Rj  ,  ...  ,  R^ 

at  points  x^  .  x^  eaually  spaced  at  a  distance  d  are  computed,  then 

each  interior  reaction  R  ,  2  S  i  %  n  -  1  is  distributed  uniformly  over  an 

interval  of  length  d  centered  at  x^  .  The  left  end  reaction  R^  is  dis¬ 
tributed  uriformly  between  x  =  0  and  x  =  0.5d  .  The  right  end  reaction  is 
distributed  similarly  between  x  =  l  -  0.5d  and  x  =  l  .  Although  this 
method  is  simple,  it  leads  to  a  small  moment  imbalance  at  the  ends.  Since  the 
load  resultant,  due  to  Rj  ,  is  moved  inward  from  the  left  end  to  a  position 
of  0.  L’r.q  ,  this  causes  a  couple  effect  which  becomes  negligible  when  enough 
supports  are  used.  A  similar  moment  imbalance  is  applicable  for  the  right 
end.  In  the  rase  where  the  loading  and  the  support  positions  are  symmetrical 
about  the  beam  center,  the  moment  imbalance,  with  respect  to  the  right  end. 


cancels  so  that  the  internal  moment  at  x  =  £  still  comes  out  to  be  zero. 

48.  A  source  listing  of  the  first  computer  program  used  to  solve  the 
three  test  examples  immediately  follows  the  numerical  results.  The  second 
program  and  numerical  results  are  detailed  and  are  given  in  Appendix  B.  Since 
both  programs  run  interactively,  the  data  input  options  are  evident  from  the 
output.  The  results  presented  were  obtained  by  use  of  an  IBM  PC/XT  computer 
employing  an  Intel-8087  floating  point  coprocessor  and  a  Winchester  disk 
drive.  All  problems  were  solved  in  less  than  two  min. 


APPENDIX  A:  PROGRAM  I — FOUNDATION  INTERACTION  BETWEEN 
A  BEAM  AND  AN  ELASTIC  HALF-PLANE 


1.  Program  I  deals  with  a  constant  depth  beam  supported  on  an  elastic 
half-plane,  in  plane  strain  or  plane  stress,  and  computes  shear,  moment, 
deflection,  and  interface  pressure.  Interaction  forces  between  the  half-plane 
and  the  beam  account  only  for  vertical  type  loads  which  occur  at  optional  and 
equally  spaced  supports  between  certain  limits  defined  according  to  the  data 


Case  I — Uniform  Distributed  Downward  Load 


2.  A  500-kip/ft  downward  pressure  load  is  applied  over  the  total  length 
of  the  bean  as  shown  in  the  following  sample  problem. 

INPUT:  Problem  title  (for  echo  check  put  S  in  column  one) 

Sample  problem  for  WES  with  constant  distributed  load  ( kips , inches ) 

SELECT  AN  OPTION:  1  =  Plane  strain,  2  =  Plane  stress 
1 

INPUT:  Young's  modulus  and  Poisson’s  ratio  for  the  half  plane 

300 .. .  2 

INPUT:  The  length,  the  depth,  and  Young's  modulus  for  the  beam 

1200. . 120. .3000. 

To  define  the  foundation  interaction  points  select:  X-min,  X-max, 
the  number  of  evenly  spaced  support  points 

and  the  width  of  the  support  in  the  direction  of  the  beam  axis 
0 . , 1200  .,41,1. 

For  the  external  beam  leading  input :  The 
number  of  concentrated  loads  and  the  number  of 
linearly  varying  ramp  loads 

0.1 

INPUT:  The  starting  magnitude,  starting  position,  end  magnitude, 
and  end  position  for  each  ramp  load 
-.A  16666666,0. . - . A16666666 , 1200 . 

Is  the  foundation  flexibility  matrix  to  be  printed?  (Y/N) 

N 

Is  the  foundation  stiffness  matrix  to  be  printed?  (Y/N) 

N 


*»*  FOUNDATION  FORCES  AND  DISPLACEMENTS  •** 


Index  Position 

Reaction 

Force 

Displacement 
in  beam 

C  relative] 

1 

. 0O000E+Q0 

1 . 58169E+01 

-6 . 83897E-14 

2 

3 . OOOOOE+Ol 

1 . 25589E+01 

-5 . 42858E-02 

3 

6 , 00000E+Q1 

1 . 14662E+01 

-1.07908E-01 

4 

9  .OOOOOE+Ol 

1 . 10904E+01 

-1 . 60276E-01 

5 

1. 20000 E+02 

1 . 10277E+01 

-2 . 1Q858E-01 

6 

1 , 50000E+02 

1 . 11146E+01 

-2. 59204E-01 

7 

1 . 80000E+02 

1 . 12720E+01 

-3.04957E-01 

8 

2 . 10000E-t-02 

1 . 14575E+01 

-3.47845E-01 

9 

2. 40000E+02 

1 . 16478E+01 

-3 . 87673E-01 

10 

2 . 70000E+02 

1 . 18295E+01 

-4 . 24309E-01 

11 

3 .  OOOOOE+02 

1 . 19957E+01 

-4 . 57679E-01 

12 

3. 30000 E+02 

1 . 21430E+01 

-4 . 87745E-01 

13 

3 . 60000E+02 

1 . 22703E+01 

-5 . 14505E-01 

14 

3 . 90000E+02 

1 . 23779E+01 

-5 . 37979E-01 

15 

4 . 20000E+02 

1 . 24669E+01 

-5 . 58200E-01 

16 

4 . 50000E+02 

1 . 25388E+01 

-5 . 75210E-01 

17 

4. 80000 E+02 

1 . 25950E+01 

-5 . 89052E-01 

18 

5 . 10000E+02 

1 . 26371E+01 

-5 . 99769E-01 

19 

5 . 40000E+02 

1 , 26661E+01 

-6.07396E-01 

20 

5. 70000E+02 

1.26832E+01 

-6 . 11961E-01 

21 

6. 00000E+02 

1 . 26888E+01 

-6. 13480E-01 

22 

6 . 30000E+02 

1.26832E+01 

-6.11961E-01 

23 

6 . 60000E+02 

1 . 26661E+01 

-6.07396E-01 

24 

6 . 90000E+02 

1 . 26371E+01 

-5.99769E-01 

25 

7. 20000E+02 

1. 25950 E+01 

-5.89052E-01 

26 

7 , 50000E+02 

1 . 25388E+01 

-5. 75210E-01 

27 

7 . 80000E+02 

1 . 24669E+01 

-5 . 58200E-01 

28 

8 , 10000E+02 

1 . 23779E+01 

-5 . 37979E-01 

29 

8 . 40000E+02 

1 . 22703E+01 

-5 . 14505E-01 

30 

8. 70000E+02 

1.21430E+01 

-4.87745E-01 

31 

9. 00000E+02 

1 . 19957E+01 

-4 . 57679E-01 

32 

9 . 30000E+02 

1. 18295E+01 

-4.24309E-01 

33 

9 . 60000E+02 

1 . 16478E+01 

-3.87673E-01 

34 

9 . 90000E+02 

1 . 14575E+01 

-3 . 47845E-01 

35 

1.02000E+03 

1 . 12720E+01 

-3 . 04957E-01 

36 

1.05000E+03 

1 . 11146E+01 

-2 . 59204E-01 

37 

1.08000E+03 

1. 10277E+01 

-2. 10858E-01 

38 

1 . 110G0E+03 

1 . 10904E+01 

-1 . 60276E-01 

39 

1.14000E+03 

1. 14662E+01 

-1.07908E-01 

40 

1 . 1 7000E+03 

1 . 25589E+01 

-5 • 42858E-02 

41 

1 . 20000E+03 

1 . 58169E+01 

.00000E+00 

Res 

idual  shear 

error  at  right  end: 

-4.08562E-14 

Residual  moment 

error  at  right  end: 

1 . 19691E-1 1 

Displacements  are  relative  to:  -1.28428E+00  at  node  #  41 


To  tabulate  foundation  pressures  for  a  given  interval:  input 
X-min,  X-max,  and  number  of  increments( input  0,0,0  to  STOP) 
0,  ,1200, ,40 


A3 


m-- 


***  SUMMARY  OF  RESULTS  FOR  INTERFACE  POINTS  *** 


Position 

Interface 

Pressure 

Shear 

Moment 

, 00000 E+00 

1.05446E+00 

.0000OE+00 

.OOOOOE+OO 

3.00000E+01 

4 . 18632E-01 

9. 59639E+00 

2. 15477E+02 

6.00000E+01 

3.82205E-01 

9. 10894E+00 

5.00155E+02 

9.00000E+01 

3. 69679E-01 

7.88721E+00 

7 . 56506E+02 

1.20000E+02 

3.67590E-01 

6.44625E+00 

9 . 71743E+02 

1 . 50000E+02 

3 . 70488E-01 

5.01741E+00 

1 . 14337E+03 

1.80000E+02 

3.75732E-01 

3 . 71071E+00 

1 . 27370E+03 

2 . 10000E+02 

3 . 81918E-01 

2. 57546E+00 

1 . 36730E+03 

2.40000E+02 

3.88259E-01 

1.62811E+00 

1 . 42964E+03 

2 . 70000E+02 

3.94317E-01 

8. 66745E-01 

1.46638E+03 

3.00000E+02 

3.99857E-01 

2. 79349E-01 

1 . 48295E+03 

3 . 30000E+02 

4 . 04766E-01 

-1 . 51311E-01 

1 . 48432E+03 

3.60000E+02 

4.09010E-01 

-4.44675E-01 

1 . 47490E+03 

3 . 90000E+02 

4 . 12597E-01 

-6 . 20570E-01 

1 . 45852E+03 

4 . 20000E+02 

4. 15565E-01 

-6 . 98137E-01 

1 . 43840E+03 

4 . 50000E+02 

4 . 17960E-01 

-6 . 95266E-01 

1 . 41723E+03 

4 . 80000E+02 

4. 19834E-01 

-6 . 28356E-01 

1 . 39717E+03 

5 . 10000E+02 

4 . 21235E-01 

-5 . 12314E-01 

1 . 37990E+03 

5.40000E+02 

4.22205E-01 

-3.60709E-01 

1 . 36670E+03 

5 . 70000E+02 

4.22774E-01 

-1 . 86027E-01 

1 . 35843E+03 

6.00000E+02 

4 . 22961E-01 

-6 . 74544E-13 

1 . 35562E+03 

6 . 30000E+02 

4 . 22774E-01 

1.86027E-01 

1. 35843E+03 

6 . 60000E+02 

4 . 22205E-01 

3.60709E-01 

1.36670E+03 

6.90000E+02 

4 . 21235E-01 

5. 12314E-01 

1 . 37990E+03 

7 . 20000E+02 

4 , 19834E-01 

6.28356E-01 

1 . 39717E+03 

7 . 50000E+02 

4 . 17960E-01 

6.95266E-01 

1 . 41723E+03 

7.80000E+02 

4 . 15565E-01 

6 . 98137E-01 

1 . 43840E+03 

8 . 10000E+02 

4 . 12597E-01 

6.20570E-01 

1 . 45852E+03 

8 . 40000E+02 

4.09010E-01 

4 . 44675E-01 

1.47490E+03 

8 . 70000E+02 

4 . 04766E-01 

1.51311E-01 

1 . 48432E+03 

9.00000E+02 

3.99857E-01 

-2, 79349E-01 

1 . 48295E+03 

9 . 30000E+02 

3.94317E-01 

-8 . 66745E-01 

1 . 46638E+03 

9 . 60000E+02 

3.88259E-01 

-1 . 62811E+00 

1 . 42964E+03 

9 . 90000E+02 

3.81918E-01 

-2 . 57546E+00 

1.36730E+03 

1 . 02000E+03 

3.75732E-01 

-3. 71071E+00 

1 . 27370E+03 

1 . 05000E+03 

3 . 70488E-01 

-5.01741E+00 

1 . 14337E+03 

1.08000E+03 

3 . 67590E-01 

-6 . 44625E+00 

9. 71743E+02 

1 . 11000E+03 

3 . 69679E-01 

-7 . 88721E+00 

7 . 56506E+02 

1 . 14000E+03 

3.82205E-01 

-9. 10894E+00 

5 . 00155E+02 

1 . 17000E+03 

4 . 18632E-01 

-9 . 59639E+00 

2. 15477E+02 

1.20000E+03 

1.05446E+00 

-1. 51545E-14 

-6 . 59668E-11 

To  tabulate  foundation  pressures  for  a  given  interval:  input 
X-min,  X-raax,  and  number  of  increments ( input  0,0,0  to  STOP) 
0,0,0 


Displacement 
[relative] 
-6. 83897E-14 
-5 , 42858E-02 
-1.07908E-01 
-1 . 60276E-01 
-2. 10858E-01 
-2 . 59204 E-01 
-3.04957E-01 
-3.47845E-01 
-3.87673E-01 
-4 . 24309E-01 
-4 . 57679E-01 
-4 . 87745E-01 
-5.14505E-01 
-5 . 37979E-01 
-5 . 58200E-01 
-5 . 75210E-01 
-5.89052E-01 
-5.99769E-01 
-6.07396E-01 
-6. 11961E-01 
-6 . 13480E-01 
-6 . 11961E-01 
-6.07396E-01 
-5. 99769E-01 
-5.89052E-01 
-5 . 75210E-01 
-5 . 58200E-01 
-5 . 37979E-01 
-5. 14505E-01 
-4 . 87745E-01 
-4 . 57679E-01 
-4 . 24309E-01 
-3.87673E-01 
-3.47845E-01 
-3.04957E-01 
-2.59204E-01 
-2. 10858E-01 
-1.60276E-01 
-1.07908E-01 
-5.42858E-02 
.00000E+00 


Case  II — Two  Concentrated  Downward  Loads  at  Beam  Ends 


3.  Two  concentrated  downward  loads  applied  at  beam  ends  appear  in  the 
following  sample  problem. 

INPUT:  Problem  title  (for  echo  check  put  S  in  column  one) 

Sample  problem  for  WES  with  concentrated  loads  at  ends(kips, inches) 

SELECT  AN  OPTION:  1  =  Plane  strain,  2  =  Plane  stress 
1 

INPUT:  Young's  modulus  and  Poisson's  ratio  for  the  half  plane 

300.. .  2 

INPUT:  The  length,  the  depth,  and  Young's  modulus  for  the  beam 

1200..  120. .3000. 

To  define  the  foundation  interaction  points  select:  X-min,  X-max, 
the  number  of  evenly  spaced  support  points 

and  the  width  of  the  support  in  the  direction  of  the  beam  axis 
0. ,1200. ,41,1. 

For  the  external  beam  loading  input :  The 
number  of  concentrated  loads  and  the  number  of 
linearly  varying  ramp  loads 

2,0 

INPUT:  The  magnitude  and  position  of  each  applied  force 
-250., 0. 

-250. ,1200. 

Is  the  foundation  flexibility  matrix  to  be  printed?  (Y/N) 
n 

Is  the  foundation  stiffness  matrix  to  be  printed?  (Y/N) 
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***  FOUNDATION  FORCES  AND  DISPLACEMENTS  **• 


Index 

Position 

Reaction 

Force 

Displacement 
in  beam 
[relative] 

1 

.00000E+00 

9 . 47707E+01 

-2 . 11667E+00 

2 

3.00000E+01 

5 . 93416E+01 

-1 . 86392E+00 

3 

6.00000E+01 

4.08453E+01 

-1 . 62026E+00 

4 

9.00000E+01 

2.88133E+01 

-1 . 39187E+00 

5 

1.20000 E+02 

2.02717E+01 

-1 . 18231E+00 

6 

1 . 50000E+02 

1 . 39668E+01 

-9 . 93318E-01 

7 

1. 80000 E+02 

9 , 23796E+00 

-8 . 25325E-01 

8 

2.10000E+02 

5.67990E+00 

-6 . 77882E-01 

9 

2. 40000E+02 

3.01605E+00 

-5 . 49948E-01 

10 

2. 70000E+02 

1.04367E+00 

-4.40118E-01 

11 

3.00000E+02 

-3 . 92848E-01 

-3 . 46791E-01 

12 

3. 30000 E+02 

-1.41 621 E+00 

-2. 68294E-01 

13 

3 , 60000E+02 

-2 . 12463E+00 

-2 . 02977E-01 

14 

3 . 90000E+02 

-2.59712E+00 

-1.492  74  E-01 

15 

4 . 20000E+02 

-2 . 89722E+00 

-1 . 05748E-01 

16 

4 . 50000E+02 

-3.07562E+00 

-7 . 11239E-02 

17 

4 . 80000E+02 

-3 . 1 7223E+00 

-4 . 43064 E-02 

18 

5 . 10000E+02 

-3 . 21772E+00 

-2. 43912E-02 

19 

5 . 40000E+02 

-3 . 23477E+00 

-1 . 0671 9E-02 

20 

5 . 70000E+02 

-3 . 23894E+00 

-2 . 64266E-03 

21 

6 . 00000E+02 

-3.23934E+00 

.00000 E+00 

22 

6 . 30000E+02 

-3 . 23894E+00 

-2 . 64266E-03 

23 

6 . 60000E+02 

-3 . 23477E+00 

-1.06719E-02 

24 

6 . 90000E+02 

-3 . 21772E+00 

-2 . 43912E-02 

25 

7 , 20000E+02 

-3. 17223E+00 

-4. 43064 E-02 

26 

7 . 50000E+02 

-3 . 07562E+00 

-7 . 11239E-02 

27 

7 . 80000E+02 

-2.89722E+00 

-1.05748E-01 

28 

8 . 10000E+02 

-2 . 59712E+00 

-1 . 49274E-01 

29 

8 . 40000E+02 

-2 . 12463E+00 

-2. 02 977 E-01 

30 

8 . 70000E+02 

-1 . 41621E+00 

-2 . 68294E-01 

31 

9. 00000 E+02 

-3 . 92848E-01 

-3 . 46791E-01 

32 

9 . 30000E+02 

1.04367E+00 

-4 . 401 18E-01 

33 

9 . 60000E+02 

3.01605E+00 

-5 . 49948E-01 

34 

9 . 90000E+02 

5 . 67990E+00 

-6 . 77882E-01 

35 

1.02000E+03 

9 . 23796E+00 

-8 . 25325E-01 

36 

1.05000E+03 

1 . 39668E+01 

-9 . 93318E-01 

37 

1.08000E+03 

2 . 0271 7E+01 

-1 . 18231E+00 

38 

1 . 11000E+03 

2 . 881 33E+01 

-1.3  9187E+00 

39 

1 . 14000E+03 

4.08453E+01 

-1 . 62026E+00 

40 

1. 17000E+03 

5 . 93416E+01 

-1 . 86392E+00 

41 

1 .20000E+03 

9 . 4  7707E+01 

-2 . 11 667E+00 

Residual  shear  error  at  right  end:  9.94760E-14 
Residual  moment  error  at  right  end:  3.99893E-11 

Displacements  are  relative  to:  -6.93861E-01  at  node  t  21 

To  tabulate  foundation  pressures  for  a  given  interval:  input 
X-rain,  X-ra<,  and  number  of  incrementsC input  0,0,0  to  STOP) 

0 .  ,  1200 .  ,  'tC 


***  SUMMARY  OF  RESULTS  FOR  INTERFACE  POINTS  *** 


Position 

.OOOOOE+OO 
3.00000E+01 
6.00000E+01 
9.00000E+01 
I . 20000E+02 
1 . 50000E+02 
1 . 80000E+02 
2. 10000E+02 
2.40000E+02 
2. 70000E+02 
3.00000E+02 
3 . 30000E+02 
3 . 60000E+02 
3 . 90000E+02 
4 . 20000E+02 
4 . 50000E+02 
4 . 8Q000E+02 
5 . 10000E+02 
5 . 4Q000E+02 
5 . 70000E+02 
6.00000E+02 
6 . 30000 E+0 2 
6 . 60000E+02 
6 . 9Q000E+02 
7 . 20000E+02 
7 . 50000E+02 
7 . 80000E+02 
8 . 10000E+Q2 
8 . 40000E+02 
8 . 70000 E+02 
9.00000E+02 
9 . 30000E+02 
9 . 60000E+02 
9 . 90000E+02 
1 . 02000E+03 
1.05000E+03 
1.08000E+03 
1.1 1000 E+0 3 
1 . 14000E+03 
1 . 17000E+03 
1 .20000E+03 


Interface 
Pressure 
6. 31804E+00 
1 . 97805E+00 
1.36L51E+00 
9. 60443E-01 
6. 75723E-01 
4. 65560E-0I 
3.07932E-01 
1 . 89330E-01 
1.00535E-01 
3. 47889E-02 
-1 . 30949E-02 
-4.72070E-02 
-7.08209E-02 
-8 . 65708E-02 
-9.65741E-02 
-1 . 02521E-01 
-1.05741E-01 
-1.07257E-01 
-1.07826E-01 
-1.07965E-01 
-1.07978E-01 
-1.07965E-01 
-1.07826E-01 
-1.07257E-01 
-1 .05741E-01 
-1.02521E-01 
-9 . 65741E-02 
-8 . 65708E-02 
-7.08209E-02 
-4 . 72070E-02 
-1 . 30949E-02 
3. 47889E-02 
1.00535E-01 
1 . 89330E-01 
3.07932E-01 
4 . 65560E-01 
6.75723E-01 
9.60443E-01 
1.36151E+00 
1 . 97805E+00 
6 . 31804E+00 


Shear 

-2.50000E+02 
-1 . 25559E+02 
-7.54650E+01 
-4.06358E+0I 
-1.60933E+01 
1 . 02599E+00 
1.26284E+01 
2.00873E+0! 
2.44353E+01 
2. 64651E+01 
2. 67905E+01 
2 . 58860E+01 
2. 41156E+01 
2. 17547E+01 
1 . 90076E+01 
1 . 60211E+01 
1 . 28972E+01 
9 . 70223E+00 
6.47599E+00 
3.23914E+00 
-1 . 12171E-12 
-3.23914E+00 
-6.47599E+00 
-9. 70223E+00 
-1.28972E+01 
-1 . 60211E+01 
-1.90076E+01 
-2. 17547E+01 
-2.4I156E+0I 
-2 . 58860E+01 
-2 . 67905E+01 
-2 . 64651E+01 
-2 . 44353E+0I 
-2.00873E+0I 
-1 . 26284E+01 
-1.02599E+00 
1 . 60933E+01 
4.06358E+01 
7 . 54650E+01 
1 . 25559E+02 
. 00000E+00 


Moment 

. OOOOOE+OO 
-5 . 14513E+03 
-8.09112E+03 
-9 . 78751E+03 
-1.06064E+04 
-1 . 08088E+04 
-1.05862E+04 
-1.00822E+04 
-9 . 40433E+03 
-8.63342E+03 
-7.82920E+03 
-7 . 03522E+03 
-6 . 28254E+03 
-5 . 59271E+03 
-4 . 980I5E+03 
-4 . 45405E+03 
-4.01991E+03 
-3 . 68075E+03 
-3. 43801E+03 
-3.29227E+03 
-3.24368E+03 
-3 . 29227E+03 
-3 . 43801E+03 
-3.68075E+03 
-4.01991E+03 
-4 . 45405E+03 
-4 . 98015E+03 
-5.59271E+03 
-6 . 28254E+03 
-7.03522E+03 
-7 . 82920E+03 
-8 . 63342E+03 
-9 . 40433E+03 
-1 . 00822E+04 
-1.05862E+04 
-1.08088E+04 
-1.06064E+04 
-9 . 78751E+03 
-8.09112E+03 
-5 . 14513E+03 
-1 . 24174E-10 


Displacement 
[relative] 
-2. 11667E+00 
-1 . 86392E+00 
-1 . 62026E+00 
-1 . 39187E+00 
-1.18231E+00 
-9 . 93318E-01 
-8 . 25325E-01 
-6 . 77882E-01 
-5.49948E-01 
-4 . 40118E-01 
-3 . 46791E-01 
-2 . 68294 E-Ol 
-2.02977E-01 
-1.49274E-01 
-1.05748E-01 
-7 . 11239E-02 
-4 . 43064E-02 
-2 . 43912E-02 
-1.06719E-02. 
-2.64266E-03 
.00000E+00 
-2.64266E-03 
-1.06719E-02 
-2.43912E-02 
-4. 43064E-02 
-7 . 11239E-02 
-1.05748E-01 
-1 . 49274E-01 
-2.02977E-01 
-2.68294E-01 
-3.46791E-01 
-4 . 40118E-01 
-5 . 49948E-01 
-6 . 77882E-01 
-8.25325E-01 
-9.93318E-01 
-1 . 18231E+00 
-1 . 39187E+00 
-1 . 62026E+00 
-1.86392E+00 
-2. 11667E+00 


To  tabulate  foundation  pressures  for  a  given  interval:  input 
X-rain,  X-max,  and  number  of  increments! input  0,0,0  to  STOP) 
0.  ,0.  ,0 


Case  III — Concentrated  End  Couples 


4.  Concentrated  end  couples  of  500  kips/ft  is  applied  at  end  of  the 
beam  as  shown  in  the  following  example. 

INPUT:  Problem  title  (for  echo  check  put  S  in  column  one) 

Sample  problem  for  WES  with  moment  on  ends  (kips , inches) 

SELECT  AN  OPTION:  1  =  Plane  strain,  2  =  Plane  stress 
1 

INPUT:  Young's  modulus  and  Poisson's  ratio  for  the  half  plane 

300.. . 2 

INPUT:  The  length,  the  depth,  and  Young's  modulus  for  the  beam 

1201..  120.. 3000. 

To  define  the  foundation  interaction  points  select:  X-min,  X-max, 
the  number  of  evenly  spaced  support  points 

and  the  width  of  the  support  in  the  direction  of  the  beam  axis 
0. ,1200. ,41,1. 

For  the  external  beam  loading  input :  The 
number  of  concentrated  loads  and  the  number  of 
linearly  varying  ramp  loads 
0,2 

INPUT:  The  starting  magnitude,  starting  position,  end  magnitude, 
and  end  position  for  each  ramp  load 
+40000. ,0. ,-40000. ,3. 

-+0000 .,1197., 40000 . , 1200 . 

Is  the  foundation  flexibility  matrix  to  be  printed?  (Y/N) 
n 


Is  the  foundation  stiffness  matrix  to  be  printed?  (Y/N) 


**♦  FOUNDATION  FORCES  AND  DISPLACEMENTS  •** 


Index 

Position 

Reaction 

Displacement 

Force 

in  beam 

1 

.00000 E+00 

-1.20259 E+02 

2.03478E+00 

2 

3.00000E+01 

-5.65107E+01 

1. 35781E+00 

3 

6.00000E+01 

-2.42488E+01 

7.97552E-01 

4 

9.00000E+01 

-5.26440E+00 

3.43477E-01 

5 

1 . 20000E+02 

6. 30331 E+00 

-1 . 67818E-02 

6 

1. 50000E+02 

1.32117E+01 

-2. 95995E-01 

7 

1 . 80000E+02 

1. 70327E+01 

-5.06591E-O1 

8 

2. 10000E+02 

1.87678E+01 

-6 . 60201E-01 

9 

2.40000E+02 

1.90987E+01 

-7 . 67416E-01 

10 

2. 70000E+02 

1.85073E+01 

-8.37669E-01 

11 

3.00000E+02 

1.73414E+01 

-8 . 79207E-01 

12 

3. 30000E+02 

1. 58539E+01 

-8 . 99128E-01 

13 

3 . 60000 E+02 

1 . 42289E+01 

-9 . 03448E-01 

14 

3. 90000 E+02 

1.25997E+01 

-8 . 97195E-01 

15 

4 . 20000E+02 

1 . 10616E+01 

-8 . 84506E-01 

16 

4 . 50000E+02 

9.68173E+00 

-8.68732E-01 

17 

4 . 80000E+02 

8. 50656E+00 

-8 . 52528E-01 

18 

5 . 10000E+02 

7.56735E+00 

-8.37945E-01 

19 

5 . 40000E+02 

6. 884 4 3 E+00 

-8.26498E-01 

20 

5. 70000E+02 

6. 47026E+00 

-8. 19225E-01 

21 

6.00000E+02 

6. 33150E+00 

-8. 16735E-01 

22 

6 . 30000E+02 

6.47026E+00 

-8. 19225E-01 

23 

6. 60000 E+02 

6.88443E+00 

-8.26498E-01 

24 

6 . 90000E+02 

7.56735E+00 

-8.37945E-01 

25 

7. 20000 E+02 

8. 50656E+00 

-8.52528E-01 

26 

7 . 50000E+02 

9.68173E+00 

-8.68732E-01 

27 

7 . 80000E+02 

1. 10616E+01 

-8.84506E-01 

28 

8.10000E+02 

1.25997E+01 

-8.97195E-01 

29 

8 . 40000E+02 

1. 42289E+01 

-9.03448E-01 

30 

8. 70000E+02 

1.58539E+01 

-8.99128E-01 

31 

9. 00000 E+02 

1. 73414E+01 

-8 . 79207E-01 

32 

9.30000E+02 

1 . 85073E+01 

-8. 37669E-01 

33 

9. 60000 E+02 

1. 90987E+01 

-7 . 67416E-01 

34 

9 . 90000E+02 

1.87678E+01 

-6.60201E-01 

35 

1.02000E+03 

1. 70327E+01 

-5.06591E-01 

36 

1.05000E+03 

1 . 32117E+01 

-2.95995E-01 

37 

1.08000E+03 

6. 30331E+00 

-1 . 67818E-02 

38 

1 , 11000E+03 

-5.26440E+00 

3. 43477E-01 

39 

1.14000E+03 

-2 . 42488E+01 

7 . 97552E-01 

40 

1 , 17000E+03 

-5.65107E+01 

1 . 35781E+00 

41 

1 . 20000E+03 

-1.20259E+02 

2.03478E+00 

Residual  shear  error  at  right  end:  .00000E+00 

Residual  moment  error  at  right  end:  -4.51337E-11 


To  tabulate  foundation  pressures  for  a  given  interval:  input 
X-rain,  X-max,  and  number  of  increment s( input  0,0,0  to  STOP) 
0., 1200., 40 
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**•  SUMMARY  OF  RESULTS  FOR  INTERFACE  POINTS  *** 


Position 

Interface 

Pressure 

Shear 

Moment 

Displacement 

.00000 E+00 

-8.01727E+00 

.00000 E+00 

.00000 E+00 

2.03478E+00 

3.00000E+01 

-1.88369E+00 

-1 . 48514E+02 

5. 70823E+04 

1.35781E+00 

6.00000E+01 

-8.08294E-01 

-1 . 88894E+02 

5. 19001E+04 

7. 97552E-01 

9.00000E+01 

-1 . 75480E-01 

-2.03651E+02 

4.59408E+04 

3.43477E-01 

1 . 20000 E+0 2 

2. 10110E-01 

-2.03131E+02 

3. 97957E+04 

-1.67818E-02 

1 . 50000E+02 

4 . 40390E-01 

-1 . 93374E+02 

3 . 38222E+04 

-2.95995E-01 

1 . 80000E+02 

5. 67755E-01 

-1 . 78252E+02 

2.82335E+04 

-5.06591E-01 

2.10000E+02 

6.25593E-01 

-1.60351E+02 

2. 31479E+04 

-6.60201E-01 

2 . 40000E+02 

6. 36623E-01 

-1 . 41418E+02 

1 . 86201E+04 

-7.67416E-01 

2. 70000E+02 

6.16911E-01 

-1.22615E+02 

1.46619E+04 

-8. 37669E-01 

3.00000E+02 

5 . 78046E-01 

-1.04691E+02 

1. 12566E+04 

-8.79207E-01 

3.30000E+02 

5.28462E-01 

-8. 80932E+01 

8. 37046E+03 

-8.99128E-01 

3 . 60000E+02 

4 . 74297E-01 

-7 . 30518E+01 

5.95938E+03 

-9.03448E-01 

3. 90000E+02 

4 . 19990E-01 

-5.96375E+01 

3. 97515E+03 

-8. 97195E-01 

4 . 20000E+02 

3.68719E-01 

-4 . 78069E+01 

2. 36925E+03 

-8.84506E-01 

4 . 50000E+02 

3 . 22724E-01 

-3 . 74352E+01 

1. 09580 E+03 

-8 . 68732E-01 

4.80000E+02 

2. 83552E-01 

-2.83411E+01 

1 . 13559E+02 

-8. 52528E-01 

5 . 10000E+02 

2 . 52245E-01 

-2.03041E+01 

-6 . 12597E+02 

-8. 37945E-01 

5.40000E+02 

2. 29481E-01 

-1.30782E+01 

-1 . 11077E+03 

-8. 26498E-01 

5. 70000 E+02 

2. 15675E-01 

-6. 40088 E+00 

-1 . 40140E+03 

-8 . 19225E-01 

6.00000E+02 

2. 11050E-01 

1 . 84741E-12 

-1.49690E+03 

-8.16735E-01 

6 . 30000E+02 

2. 15675E-01 

6. 40088E+00 

-1.40140E+03 

-8. 19225E-01 

6. 60000E+02 

2.29481E-01 

1.30782E+01 

-1 . 11077E+03 

-8.26498E-01 

6 . 90000E+02 

2.52245E-01 

2.03041E+01 

-6. 12597E+02 

-8. 37945E-01 

7 . 20000E+02 

2.83552E-01 

2.83411E+01 

1.13559E+02 

-8. 52528E-01 

7 . 50000E+02 

3 . 22724E-01 

3. 74352E+01 

1.09580E+03 

-8 . 68732E-01 

7 . 80000E+02 

3.68719E-01 

4 . 78069E+01 

2. 36925E+03 

-8.84506E-01 

8 . 10000E+02 

4 . 19990E-01 

5 . 96375E+01 

3. 97515E+03 

-8.97195E-01 

8.40000E+02 

4. 74297E-01 

7. 30518E+01 

5 . 95938E+03 

-9.03448E-01 

8 . 70000E+02 

5.28462E-01 

8 . 80932E+01 

8. 37046E+03 

-8. 99128E-01 

9.00000E+02 

5 . 78046E-01 

1.04691E+02 

1 . 12566E+04 

-8. 79207E-01 

9 . 30000E+02 

6 . 16911E-01 

1 . 22615E+02 

1.46619E+04 

-8 . 37669E-01 

9 . 60000E+02 

6. 36623E-01 

1 . 41418E+02 

1 . 86201E+04 

-7.67416E-01 

9 . 90000E+02 

6. 25593E-01 

1 . 60351E+02 

2. 31479E+04 

-6.60201E-01 

1.02000E+03 

5.67755E-01 

1.78252E+02 

2.82335E+04 

-5.06591E-01 

1.05000E+03 

4 . 40390E-01 

1.93374E+02 

3.38222E+04 

-2.95995E-01 

1.08000E+03 

2. 10110E-01 

2.03131E+02 

3. 97957E+04 

-1.67818E-02 

1 . 11000E+03 

-1 . 75480E-01 

2.03651E+02 

4 . 59408E+04 

3.43477E-01 

1 . 14000E+03 

-8.08294E-01 

1.88894E+02 

5. L9001E+04 

7.97552E-01 

1 . 17000E+03 

-1 . 88369E+00 

1 . 48514E+02 

5. 70823E+04 

1 . 35781E+00 

1.20000E+03 

-8.01727E+00 

3.60766E-12 

-9 . 96564E-10 

2.03478E+00 

To  tabulate  foundation  pressures  for  a  given  interval:  input 
X-min,  X-raax,  and  number  of  increraents( input  0,0,0  to  STOP) 
0. ,0.  ,0 


5.  The  computer  codes  for  Program  I  sample  problems  are  shown  on  the 
following  pages. 

1:$$$$$$$$  MAIN 
2:C 

3:C  Written  for:  The  U.S.  Corps  of  Engineers 

4:C  By:  H.B.  Wilson  and  L.H.  Turcotte 

5:C  University  of  Alabama  at  Tuscaloosa 

6:C  Date:  September  1983 

7 :  C 

8:  IMPLICIT  REAL*8 

9:  CHARACTER* 1 

10:  LOGICAL 

11:  COMMON  /  BEAM  / 

12:  & 

13:  & 

14:  COMMON  /  QTOTAL  / 

15:  COMMON  /  HAFPLA  / 

16:  COMMON  /  EQNS  / 

17:  COMMON  /  EQNS2  / 

18:  COMMON  /  EQNS3  / 

19:  DATA 

20:  DATA 

21 :  C 

22:  NFMAX  =  53 

23:  ECHO  =  .FALSE. 

24:  IFRST  =  1 

25 :  C 

26:  WRITE  (IPRT, 10) 

27:  10  FORMAT ( // , 5X , ' *** 

28:  &  /,5X, ’ *** 

29:  WRITE  (IPRT, 20) 

30:  20  FORMAT ( / , '  INPUT: 

31:  &  ’one)') 

32:  READ  (ICRT.30)  (TITL(I),I  =  1,80) 

33:  30  FORMAT (80A1) 

34:  IF  (TITL(l)  .EQ.  '$')  ECHO  =  .TRUE. 

35:  IF  (ECHO)  WRITE  (IPRT,40)  (TITL(I),I  =  2,80) 

36:  40  FORMAT (IX, 80A1) 

37:  WRITE  (IPRT.50) 

38:  50  FORMATt/,'  SELECT  AN  OPTION:  1  =  P^ane  strain,  2  =  Plane  stress') 

39:  READ  (ICRT,*)  ISTATE 

40:  IF  (ECHO)  WRITE  (IPRT,*)  ISTATE 

41:  WRITE  (IPRT, 60) 

42:  60  FORMAT (/,'  INPUT:  Young"  s  modulus  and  Poisson' 's  ratio  for  ', 

43:  &  'the  half  plane') 

44:  READ  (ICRT,*)  EHAPL.POIHPL 

45:  IF  (ECHO)  WRITE  (IPRT,*)  EHAPL.POIHPL 

46:C. . .Define  elastic  constants 
47:  TWOGHP  *  EHAPL  /  (1.D0  +  POIHPL) 

48 :C... Plane  strain 

49:  IF  (ISTATE  .EQ.  1)  CAPA  =  3. DO  -  4. DO  *  POIHPL 

50:C... Plane  stress 

51:  IF  (ISTATE  .EQ.  2)  CAPA  =  (3. DO  -  POIHPL)  /  (1.D0  +  POIHPL) 


(A-H.O-Z) 

TITL(8Q) ,ANS 
ECHO 

BLEN.XFIX.NF, F(53 ) ,XF( 53 ) ,NP,P(53) , 
XP(53) ,NR,WL(53) ,XL(53) ,WR(53) ,XR(53) , 
EIB 

XT( 53), IFRST, D.DH 

EHAPL , POIHPL , TWOGHP , CAPA 

NFMAX,AHP(53,53) ,ABM( 53,53 ),B( 53) 

ATOTL( 53 , 53 ) , IPIV0T( 53 ) , SCAL( 53 ) , ISTATE 
AHPINV ( 5  3 , 53 ) , DBM ( 5  3 ) , SMALL , NODEI 
PI/  3. 14 1592653 589 79 3D0  / 

IRD,  IPRT  /  0,  0  / 


FOUNDATION  INTERACTION  BETWEEN  A  BEAM  ***’ 

AND  AN  ELASTIC  HALF  PUNE  ***',/) 

Problem  title  (for  echo  check  put  $  in  column  ’, 


All 


52:  SHRMOD  =  EHAPL  •  .5DO  /  (l.DO  +  POIHPL) 

53:  WRITE  (IPRT, 70) 

54:  70  FORMAT (/,'  INPUT:  The  length,  the  depth,  and  Young' 's  modulus  ', 

55:  &  '  for  the  beam') 

56:  READ  (ICRT,*)  BLEN , BDEP , EBEAM 

57:  IF  (ECHO)  WRITE  (IPRT,*)  BLEN , BDEP , EBEAM 

58:  EIB  =  EBEAM  •  (BDEP  **  3)  /  12. DO 

59:  XFIX  =  1.1D0  *  BLEN 

60:  WRITE  (IPRT, 80) 

61:  80  F0RMAT(/,'  To  define  the  foundation  interaction  points  select:  ', 

62:  &  'X-min,  X-max,',/,'  the  number  of  evenly  spaced  support', 

63:  &  '  joints',/,'  and  the  width  of  the  support  in  the  ', 

64:  &  'direction  of  the  beam  axis') 

65:  READ  (ICRT,*)  XMIN , XMAX, NF , WDTH 

66:  IF  (ECHO)  WRITE  (IPRT,*)  XMIN, XMAX, NF, WDTH 

67:  IF  (NF  .LE.  NFMAX-2)  GO  TO  100 

68:  WRITE  (IPRT, 90)  NFMAX-2 

69:  90  F0RMAT(/,'  Maximum  number  of  support  points  =  ',13, 

70 :  &  '  PROGRAM  TERMINATED ' ) 

71:  STOP  '  ' 

72:  100  DF  =  (XMAX  -  XMIN)  /  (NF  -  1) 

73:  HAFWID  =  0.5DO  *  WDTH 

74:  ALP  =  (l.DO  +  CAPA)  /  (8. DO  •  HAFWID  *  SHRMOD  *  PI) 

75:  DO  110  I  =  1 , NF 

76:  110  XF( I )  =  XMIN  +  DF  *  (I  -  1) 

77:  IF  (BLEN  .EQ.  0.0D0)  GO  TO  180 

78:  WRITE  (IPRT, 120) 

79:  120  F0RMAT(/, '  For  the  external  beam  loading  input:  The  ' ,/, 

80:  &  '  number  of  concentrated  loads  and  the  number  of  ' ,/, 

81:  &  '  linearly  varying  ramp  loads  ’) 

82:  READ  (ICRT,*)  NP.NR 

83:  IF  (ECHO)  WRITE  (IPRT,*)  NP,NR 

84:  IF  (NP  .EQ.  0)  GO  TO  150 

85:  WRITE  (IPRT, 130) 

86:  130  FORMAT (/,'  INPUT:  The  magnitude  and  position  of  each  applied  ', 

87 :  &  ' force ' ) 

88:  DO  140  I  =  1,NP 

89:  READ  (ICRT,*)  P(I),XP(I) 

90:  IF  (ECHO)  WRITE  (IPRT,*)  P(I),XP(I) 

91:  140  CONTINUE 

92:  150  IF  (NR  .EQ.  0)  GO  TO  180 
93:  WRITE  (IPRT, 160) 

94:  160  F0RMAT(/,'  INPUT:  The  starting  magnitude-,  starting  ’, 

95:  &  'position,  end  magnitude,',/,'  and  end  position  for  ', 

96:  &  'each  ramp  load') 

97:  DO  170  I  =  1,NR 

98:  READ  (ICRT,*)  WL(I) ,XL(I) ,WR(I) ,XR(I> 

99:  IF  (ECHO)  WRITE  (IPRT,*)  WL(I) ,XL(I) ,WR(I) ,XRvI) 

100:  170  CONTINUE 

101 :C... Form  equations  and  solve  for  the  foundation  equations 
102:  180  CALL  S0LEQN( HAFWID, ALP ) 

103:  WRITE  (IPRT, 190) 

104:  190  F0RMAT(/,'  Is  the  foundation  flexibility  matrix  to  be', 

105:  &  '  printed?  (Y/N)') 


A12 


106:  READ  (ICRT.30)  ANS 

107:  IF  (ECHO)  WRITE  (IPRT,*)  ANS 

108:  IF  (ANS  .EQ.  'N'  .OR.  ANS  .EQ.  ’n’)  GO  TO  230 

109:  WRITE  (IPRT, 200) 

110:  200  FORMAT (/,’  The  Foundation  Flexibility  Matrix  is:’,/) 

HI:  DO  210  I  =  1  ,NF 

112;  210  WRITE  (IPRT.220)  (AHP(I,J),  J  =  l.NF) 

113:  220  FORMAT ( IX, 6 (1PE 12. 4) ,/,(13X,5(lPE12.4))) 

114:  230  WRITE  (IPRT.240) 

115:  240  FORMAT(/,'  Is  the  foundation  stiffness  matrix  to  be  printed?  ' 
116:  &  '(Y/N)') 

117:  READ  (ICRT.30)  ANS 

118:  IF  (ECHO)  WRITE  (IPRT,*)  ANS 

119:  IF  (ANS  .EQ.  'N'  .OR.  ANS  .EQ.  'n')  GO  TO  290 

120:  CALL  INVERT ( AHP , NF , NFMAX , IPI VOT , SCAL , IFLAG , AHPINV ) 

121:  IF  (IFLAG  .NE.  2)  GO  TO  260 

122:  WRITE  (IPRT, 250) 

123:  250  FORMAT(/,'  The  Foundation  Flexibility  Matrix  is  Singular 
124:  &  /,'  EXECUTION  TERMINATED* ,/) 

125:  STOP  '  ’ 

126:  260  WRITE  (IPRT.270) 

127:  270  FORMAT(/,'  The  Foundation  Stiffness  Matrix  is:',/) 

128:  DO  280  I  =  l.NF 

129:  280  WRITE  (IPRT.220)  (AHPINV(I, J) , J  =  l.NF) 

130:  290  CONTINUE 

131:  WRITE  (IPRT, 300) 

132:  300  FORMAT (/,3X, '  *»*  FOUNDATION  FORCES  AND  DISPLACEMENTS  •••* , 
133:  &  //,'  Index' ,2X, 'Position' ,8X, 'Reaction' ,6X, 

134:  &  'Displacement' ,/,24X, ' Force’ ,11X, ' in  beam’) 

135:  IF  (SMALL  .NE.  0.D0)  WRITE ( IPRT, 310) 

136.  310  F0RMAT(39X,’ [relative] ') 

137:  SUMF  =  Q(BLEN,1) 

138:  SUMM  =  Q( BLEN ,2) 

139:  DO  320  1*1, NF 

140:  SUMF  =  SUMF  -  F(I) 

141:  SUMM  =  SUMM  -  F(I)  *  (BLEN-XF(I) ) 

142:  320  WRITE( IPRT, 330)  I,XF(I) ,-F(I) ,DBM(I) 

143:  330  FORMAT( 14 , 3( 1PE13. 5 , 3X) ) 

144:  WRITE ( IPRT, 340)  SUMF, SUMM 

145:  340  FORMAT(/,'  Residual  shear  error  at  right  end:  ’.1PE13.5, 

146:  &  /,'  Residual  moment  error  at  right  end:  ',1PE13.5) 

147:  IF  (SMALL  .NE.  O.DO)  WRITE( IPRT, 350 )  SMALL, NODEI 

148:  350  FORNAT(/,'  Displacements  are  relative  to<  ',1PE13.5, 

149:  &  '  at  node  #  ',  13) 

150 : C. . .Compute  shear  and  moment  at  interface  points  and  output  results 
151:  360  WRITE (IPRT, 3 70) 

152:  370  F0RMAT(//,'  To  tabulate  foundation  pressures  for  a  given  ', 

153:  &  'interval:  input',/,'  X-min,  X-max,  and  number  ', 

154:  &  'of  increments ( input  0,0,0  to  STOP)') 

155:  READdRD,*)  XXL, XXR, NSEC 

156:  IF  (XXL  .EQ,  XXR)  STOP  '  ' 

157:  DX  =  (XXR- XXL)  /  DBLE(NSEG) 

158:  WRITE (IPRT, 380) 

159:  380  FORMAT (/, 15X, ' •**  SUMMARY  OF  RESULTS  FOR  INTERFACE  POINTS  ’, 


&  //, 

&  4X, 'Position' ,7X, 'Interface' ,8X, 'Shear' ,9X, 'Moment 

&  6X, ' Displacement ' ) 

IF  (SMALL  .NE.  0.0D0)  WRITE< IPRT , 390 ) 

390  F0RMATQ9X,  'Pressure'  ,36X, '  [relative] ' ) 

IF  (SMALL  .EQ.  0.0D0)  WRITE (IPRT, 400) 

400  FORMAT* 19X, ' Pressure' ) 

DO  410  I  =  O.NSEG 
XX  =  XXL  +  I  •  DX 
QT1  =  QTOTL(XX.l) 

QT2  =  QT0TL(XX, 2) 

QT4  =  QT0TL(XX,4) 

FPR  =  FPRES(XX) 

410  WRITE ( IPRT, 420)  XX,FPR,QT1,QT2,QT4 
420  F0RMAT(1X,5(1PE13.5,2X) ) 

GO  TO  360 
END 

C$$$$$$$$$$  FS 

REAL*8  FUNCTION  FS(X,XBEGIN,NPOWR) 

C. . .Singularity  function 

IMPLICIT  REAL*8  (A-H.O-Z) 

FS  =  O.ODO 

D  =  X  -  XBEGIN 

IF  (D  .LT.  O.ODO)  RETURN 

IF  ( NPOWR  .GT.  0)  GO  TO  10 

FS  =  1 . DO 

RETURN 

10  FS  =  D  ••  NPOWR 
RETURN 
END 

CSSSSSSSSSS  RAMP 

REAL*8  FUNCTION  RAMP(X,XL,WL,XR,WR,EI,ID) 

C...Load  and  deflection  quantities  for  ramp  load  on  a  beam 
IMPLICIT  REAL*8  (A-H.O-Z) 

RAMP  =  O.ODO 
DL  =  X  -  XL 

IF  (DL  .LT.  O.ODO)  RETURN 
DR  =  X  -  XR 

S  =  (WR  -  WL)  /  (XR  -  XL) 

C 

GO  TO  (10,20,30,40,50), ID+1 
C 

C...Load  per  unit  length 

10  RAMP  =  WL  •  FS(X,XL,0)  -  WR  •  FS(X,XR,0) 

IF  (S  .EQ.  O.ODO)  RETURN 

RAMP  =  RAMP  +  S  •  (FS(X,XL,1)  -  FS(X,XR,1)) 

RETURN 
C . . . Shear 

20  RAMP  =  WL  *  FS(X,XL, 1 )  -  WR  *  FS(X,XR,1) 

IF  (S  .EQ.  O.ODO)  RETURN 

RAMP  =  RAMP  +  0.5D0  *  S  *  (FS(X,XL,2)  -  FS(X,XR,2)) 

RETURN 
C. . .Moment 

30  RAMP  =  0.5D0  *  (WL  *  FS(X,XL,2)  -  WR  *  FS(X,XR,2)) 


214 :  IF  (S  .EQ.  O.ODO)  RETURN 

215:  RAMP  =  RAMP  +  S  /  6.0D0  *  (FS(X,XL,3)  -  FS(X,XR,3)) 

216:  RETURN 

217 : C. . . Slope 

218:  40  RAMP  =  (WL  *  FS(X,XL,3)  -  HR  •  FS(X,XR,3))  /  6.0D0 

219:  IF  (S  .EQ.  O.ODO)  GO  TO  60 

220:  RAMP  =  RAMP  +  S  /  24.0D0  *  (FS(X,XL,4)  -  FS(X,XR,4)) 

221:  GO  TO  60 

222 :C. . . Deflect ion 

223:  50  RAMP  =  (WL  *  FS(X,XL,4)  -  WR  *  FS(X,XR,4))  /  24.0D0 

224:  IF  (S  .EQ.  O.ODO)  GO  TO  60 

225:  RAMP  =  RAMP  +  S  /  120. ODO  *  (FS(X,XL,5)  -  FS(X,XR,5)) 

226:  60  RAMP  =  RAMP  /  El 

227:  RETURN 

228:  END 

229:0$$$$$$$$$$  FUNIT 

230:  REAL*8  FUNCTION  FUNIT ( X , POSITN , El , ID ) 

231:C...Load  and  deflection  quantities  for  concentrated  load  on  a  beam 
232:  IMPLICIT  REAL*8  (A-H.O-Z) 

233:  FUNIT  =  O.ODO 

234:  IF  (ID  .EQ.  0)  RETURN 

235:  D  =  X  -  POSITN 

236:  IF  (D  .LT.  O.ODO)  RETURN 

237:  GOTO  (10,20,30,40)  , ID 

238 :C. . .Shear 

239:  10  FUNIT  =  FS(X, POSITN, 0) 

240:  RETURN 

241 : C. . .Moment 

242:  20  FUNIT  =  FS(X, POSITN, 1) 

243:  RETURN 

244 :C. . . Slope 

245:  30  FUNIT  =  FS(X , POSITN , 2 )  *  0.5D0  /  El 

246:  RETURN 

247 :C. . .Deflection 

248:  40  FUNIT  =  FS(X, POSITN, 3)  /  (6. ODO  *  El) 

249:  RETURN 

250 :  END 

251 :C$$$$$$$$$$  Q 

252:  REAL*8  FUNCTION  Q(X,ID) 

253 :C. .. Loading  and  deflection  contribution  of  external  concentrated 

254 :C... and  ramp  loads 

255:  IMPLICIT  REAL*8  (A-H.O-Z) 

256:  COMMON  /  BEAM  /  BLEN ,XFIX, NF, F( 53 ) ,XF(53 ) ,NP , P( 53 ) , 

257:  &  XP(53),NR,WL(53),XL(53),WR(53),XR(53), 

258:  &  EIB 

259:  Q  =  O.ODO 

260:  IF  (NP  .EQ.  0)  GO  TO  20 

261:  DO  10  J  =  1 ,NP 

262:  10  Q  =  Q  +  P(J)  •  FUNIT (X,XP(J) , EIB, ID) 

263:  20  IF  (NR  .EQ.  0)  RETURN 

264:  DO  30  J  =  1,NR 

265:  30  Q  =  Q  +  RAMP(X,XL( J) ,WL(J) ,XR( J) ,WR( J) , EIB, ID) 

266:  RETURN 

267:  END 


268 : CS$$$$$$$$$  QTOTL 

269:  REAL* 8  FUNCTION  QTOTL (X, ID) 

270:C. . .Total  loading  and  deflection  quantities 
271:  IMPLICIT  REAL*8  (A-H.O-Z) 

272:  COMMON  /  BEAM  /  BLEN,XFIX,NF,F(53) ,XF(53) ,NP,P(53) , 

273:  &  XP(53) ,NR,WL(53) ,XL(53) ,WR(53) ,XR(53) , 

274:  &  EIB 

275:  COMMON  /  QTOTAL  /  XT ( 53 ) , IFRST , D , DH 

276:  COMMON  /  EQNS3  /  AHPINV(53 , 53) ,DBM( 53) , SMALL, NODEI 

277 :C. .. Initialize  array  XT  on  first  entry 

278:  IF  (IFRST  .EQ.  0)  GO  TO  20 

279:  IFRST  =  0 

280:  N1  =  NF  -  1 

281:  XT(1)  =  XF(1) 

282:  XT(NF+1 )  =  XF(NF) 

283:  D  =  (XF( NF)-XF( 1 ) )/DBLE(Nl ) 

284 :  DH  *  0.5D0  *  D 

285:  DO  10  I  =  2,  NF 

286:  10  XT( I )  =  XF(I)  -  DH 

287 : C . . . 

288:  20  IF  (ID  .GT.  2)  GO  TO  50 

289:  XLT  =  XT(1) 

290:  XRT  *  XT(2) 

291:  WT  =  2.0D0  *  F(l) 

292:  QT  =  RAMP(X,XLT,WT,XRT,WT,EIB,ID) 

293:  DO  30  I  =  2,  N1 

294:  XLT  =  XT(I) 

295:  IF  (XLT  .GE.  X)  GO  TO  40 

296:  XRT  *  XT(I+1) 

297:  WT  =  F(I) 

298:  30  QT  =  QT  +  RAMP(X,XLT,WTtXRT,WT,EIB, ID) 

299:  XLT  =  XT(NF) 

300:  IF  (X  .LE.  XLT)  GO  TO  40 

301:  XRT  =  XT(NF+1 ) 

302:  WT  =  2.0D0  *  F(NF) 

303:  QT  =  QT  +  RAMP(X,XLT,WT,XRT,WT,EIB,ID) 

304:  40  QTOTL  =  Q(X,ID)  -  QT  /  D 

305 :  RETURN 

306 : C. .. Interpolate  lineraly  for  slope  or  deflection 
307:  50  ID2  =  ID  -  2 

308:  QTOTL  =  FLNTRP ( X , XF , DBM , NF , ID2 ) 

309 :  RETURN 

310:  END 

311 :C$$$$$$$$$$  FPRES 

312:  REAL*8  FUNCTION  FPRES(X) 

313 :C. .. Function  to  determine  piecewise  constant  foundation  pressure 
314:  IMPLICIT  REAL*8  (A-H.O-Z) 

315:  COMMON  /  QTOTAL  /  XT( 53 ), IFRST , D,DH 

316:  COMMON  /  BEAM  /  BLEN,XFIX,NF,F(53) ,XF(53) ,NP,P(53) , 

317:  &  XP(53) ,NR,WL(53) ,XL(53) ,WR(53) ,XR(53) , 

318:  &  EIB 

319:  IF  (X  .GT.  XT(2))  GO  TO  10 

320:C. . .Pressure  in  short  left  segment 

321:  FPRES  =  2. DO  *  F(l) 


•ft; 


n*%,  / 


A16 


322:  GO  TO  50 

323:  10  IF  (X  .LT.  XT(NF) )  GO  TO  20 

324 : C. . .Pressure  in  right  short  segment 
325:  FPRES  =  2. DO  *  F(NF) 

326:  GO  TO  50 

327 : C. . .Determine  segment  containing  X 
328:  20  DO  30  I  =  3,NF 

329:  ISV  =  I 

330:  IF  ( XT ( I )  .GE.  X)  GO  TO  40 

331:  30  CONTINUE 

332:  40  FPRES  =  F(ISV-l) 

333:  50  FPRES  =  -FPRES  /  D 

334:  RETURN 

335:  END 

336 : C$$$$$$$$S$  FACTR 

337:  SUBROUTINE  FACTR  (A, N, NMAX, IPIVOT, SCAL, IFLAG) 

338:C. . .Perform  triangular  factorization  of  matrix  a 

339 : C. . .using  scaled  row  pivoting 

340:C...  IFLAG  =  1  means  normal  return 

341 :C...  =  2  means  matrix  is  singular 

342:  IMPLICIT  REAL*8  (A-H,0-Z) 

343:  DIMENSION  A(NMAX,N) ,  IPIVOT(N),  SCAL(N) 

344:  IFLAG  =  1 

345:C. . .Initialize  IPIVOT  and  SCAL 

346:  DO  20  I  =  1,N 

347:  IPIVOT(I)  =  I 

348:  ROWMAX  =  0.0D0 

349:  DO  10  J  *  1,N 

350:  10  ROWMAX  =  DMAX1 (ROWMAX, DABS (A( I ,J) ) ) 

351:  IF  ( ROWMAX. EQ.0.0D0)  GO  TO  50 

352:  20  SCAL(I)  =  ROWMAX 

353:C. . .Perform  Gauss  reduction 

354:  NM1  =  N-l 

355:  IF  (NM1.EQ.0)  RETURN 

356:  DO  40  K  =  l.NMl 

357:  J  =  K 

358:  KP1  =  K+l 

359:  IP  =  IPIVOT(K) 

360:  COLMAX  =  DABS(A(IP,K) ) /SCAL (IP) 

361:  DO  30  I  =  KPl.N 

362:  IP  =  IPIVOT ( I ) 

363:  AWIKOV  =  DABS( A( IP, K) )/SCAL( IP ) 

364:  IF  ( AWIKOV. LE. COLMAX)  GO  TO  30 

365:  COLMAX  =  AWIKOV 

366:  J  =  I 

367:  30  CONTINUE 

368:  IF  ( COLMAX. EQ.O.ODO)  GO  TO  50 

369:  IPK  =  IPIVOT(J) 

370:  IPIVOT(J)  =  IPIVOT(K) 

371:  IPIVOT ( K )  =  IPK 

372:  DO  40  I  =  KPl.N 

373:  IP  =  IPIVOT ( I ) 

374:  A(IP,K)  =  A(IP,K)/A(IPK,K) 

375:  RATIO  =  -A(IP,K) 


376:  DO  40  J  *  KPl.N 

377:  40  A(IP,J)  *  RATIO*A(IPK, J)+A(IP, J) 

378:  IF  (A(IP.N).EQ. O.ODO)  GO  TO  50 

379 :  RETURN 

380:  50  IFLAG  -  2 

381:  RETURN 

382 :  END 

383 : CSSSSSSSSSS  BAKSUB 

384:  SUBROUTINE  BAKSUB  (A,N,NMAX,B,IPIVOT,X) 

385:C. . .Solve  simultaneous  equations  AX  -  B  where  matrix  A  has 
386 :C... been  subjected  to  factorization  by  SUBROUTINE  FACTR 
387:  IMPLICIT  REAL*8  (A-H.O-Z) 

388:  DIMENSION  A(NMAX,N) ,  B(N) ,  IPIVOT(N),  X(N) 

389:  IF  (N.GT.l)  GO  TO  10 

390:  X(l)  =  B( 1 )/A( 1,1) 

391:  RETURN 

392:  10  IP  =  IPIVOT(l) 

393:  X(l)  =  B( IP) 

394:  DO  30  K  *  2,N 

395:  IP  =  IPIVOT(K) 

396:  KM1  «  K-l 

397:  SUM  =  O.ODO 

398;  DO  20  J  =  l.KMl 

399:  20  SUM  *  A(IP, J)*X( J)+SUM 

400:  30  X(K)  *  B(IP)-SUM 

401:  X(N)  =  X(N)/A(IP,N) 

402:  K  =  N 

403:  DO  50  NP1MK  =  2,N 

404:  KP1  »  K 

405:  K  *  K-l 

406:  IP  =  IPIVOT(K) 

407:  SUM  =  O.ODO 

408:  DO  40  J  =  KPl.N 

409:  40  SUM  =  A(IP, J)*X( J)+SUM 

410:  50  X(K)  =  (X(K)-SUM)/A(IP,K) 

411:  RETURN 

412:  END 

413 :C$$$$$$$$$$  INVERT 

414:  SUBROUTINE  INVERT(A,N,NMAX,IPIVOT,SCAL,IFLAGtAINV) 

415 :C. .. IFLAG  indicates  singular  matrix 
416:  IMPLICIT  REAL*8  (A-H.O-Z) 

417:  DIMENSION  A(NMAX, 1 ) , AINV(NMAX, 1 ) , IPIV0TC1 ) , SCALC1 ) 

418:  CALL  FACTR( A.N.NMAX.IPIVOT, SCAL.IFLAG)  . 

419:  IF  (IFLAG  .EQ.  2)  RETURN 

420:  DO  10  J  =  l.N 

421:  10  SCAL(J)  *  O.ODO 

422:  DO  20  I  =  l.N 

423:  SCAL(I)  =  l.ODO 

424:  CALL  BAKSUB ( A.N.NMAX, SCAL, IPIVOT , AINV(1 , I) ) 

425:  SCAL(I)  =  O.ODO 

426:  20  CONTINUE 

427:  RETURN 

428:  END 

429: CSSSSSSSSSS  SOLEQN 


430:  SUBROUTINE  SOLEQN ( HAFW ID , ALP ) 

431:C...This  subroutine  forms  the  flexibility  matrices  for  a  beam 
432:C...and  for  a  half  plane  subjected  to  several  concentrated 
433:C. . .loads.  These  matrices  are  then  combined  and  the  resulting 
434 : C. . .simultaneous  equations  are  solved  for  the  interaction  forces. 
435:  IMPLICIT  REAL* 8  (A-H.O-Z) 

436:  COMMON  /  BEAM  /  BLEN,XFIX,NF,F(53) ,XF(53) ,NP,P(53) , 

437:  &  XP(53) ,NR,WL(53) ,XL(53) ,WR(53) ,XR(53) , 

438:  &  '  EIB 

439:  COMMON  /  HAFPLA  /  EHAPL.POIHPL.TWOGHP.CAPA 

440:  COMMON  /  EQNS  /  NFMAX, AHP ( 53 , 53 ) ,ABM(53,53) ,B(53) 

441:  COMMON  /  EQNS2  /  AT0TL( 53 , 53 ) , IPIV0T(53 ) ,SCAL( 53 ) , ISTATE 

442:  COMMON  /  EQNS3  /  AHPINV(53, 53) ,DBM(53) .SMALL, NODEI 

443:  DATA  IPRT/  0  / 

444:  NF1  =  NF  +  1 

445:  NF2  =  NF  +  2 

446:C. . .Define  some  scaling  factors 
447:  SCAL1  =  BLEN  **  3  /  EIB 

448:  SCAL2  =  BLEN  **  2  /  EIB 

449 :C... Form  influence  for  half  plane 
450:  DO  10  I  =  1,NF 

451:  DO  10  J  =  1 ,NF 

452:  10  AHP(I.J)  =  FORIGN (XF(I)-XF(J), HAFWID , XFIX , ALP ) 

453 :C... Form  beam  influence  if  requested 
454:  IF  (BLEN  .EQ.  O.ODO)  RETURN 

455:  DO  20  I  =  1,NF 

456:  DO  20  J  =  1,NF 

457:  20  ABM(I.J)  =  FS(XF(I) ,XF(J) ,3)  /  (6.0D0  *  EIB) 

458 :C... Form  right  side  vector  for  applied  beam  loads 
459:  DO  30  I  *  1,NF 

460:  30  B(I)=Q(XF(I) ,4) 

461 :C... Form  combined  influence 

462:  DO  40  I  =  1,NF 

463:  DO  40  J  =  1,NF 

464:  40  ATOTL(I.J)  =  AHP(I.J)  +  ABM(I.J) 

465:C. . .Define  additional  static  equilibrium  equations 

466:  DO  50  J  =  1,NF 

467:  ATOTL(J.NFl)  =  -SCAL1 

468:  ATOTL ( J , NF2 )  =  -  XF(J)  *  SCAL2 

469:  ATOTL(NFl.J)  =  1.0D0 

470:  50  ATOTL (NF2.J)  =  BLEN  -  XF(J) 

471:  ATOTL (NF1.NF1)  =  O.ODO 

472:  ATOTL (NF1.NF2)  =  O.ODO 

473:  ATOTL (NF2.NF1)  =  O.ODO 

474:  AT0TL(NF2,NF2)  =  O.ODO 

475:  B(NF1)  =  Q(BLEN.l) 

476:  B(NF2)  =  Q(BLEN,2) 

477 : C. .. Solve  for  the  interaction  forces 

478:  CALL  FACTR( ATOTL, NF2, NFMAX, IPIVOT, SCAL.IFLAG) 

479:  IF  ( I FLAG  .NE.  2)  GO  TO  70 

480:  WRITE  (IPRT, 60) 

481:  60  F0RMAT(/,’  COMBINED  STIFFNESS  MATRIX  IS  SINGULAR’,/, 

482:  &  ’  EXECUTION  IS  TERMINATED’) 

483:  STOP  ’  ’ 


484 :C. . .Compute  the  forces  at  interaction  points 
485:  70  CALL  BAKSUBC ATOTL ,  NF2 , NFMAX , B , IPI VOT , F ) 

486 : C. . .Compute  left  end  slope  and  displacement 
487:  F(NF2)  =  F(NF2)  *  SCAL2 

488:  F(NF1)  =  F(NF1)  *  SCAL1 

489 : C. . .Compute  the  displacements  at  interaction  points 
490:  DO  90  I  =  1,NF 

491:  DBMI  =  B(I)  +  F(NF1)  +  F(NF2)  *  XF(I) 

492:  DO  80  J  =  1,NF 

493:  80  DBMI  =  DBMI  -  ABM(I.J)  *  F(J) 

494:  90  DBM(I)  =  DBMI 

495:C. . .Find  smallest  and  largest  displacement  in  beam 

496:  BIC  =  +1.D20 

497:  SMALL  =  -1.D20 

498:  DO  100  I  =  1,NF 

499:  IF (  DBM(I)  .CT.  SMALL  )  THEN 

500:  SMALL  =  DBM(I) 

501:  NODEI  =  I 

502:  ENDIF 

503:  IF (  DBM(I)  .LT.  BIG  )  THEN 

504:  BIG  =  DBM(I) 

505:  ENDIF 

506 :  100  CONTINUE 

507 :C... Make  all  displacements  relative  to  SMALL=0  if  displacements  do 
508:C...not  change  sign 

509:  IF(  BIG*SMALL  .LT.  0.D0  )  SMALL  *  0.D0 

510:C. . .Shift  all  displacements  by  value  of  SMALL 

511:  DO  110  I  =1 , NF 

512;  110  DBM(I)  =  DBM(I)  -  SMALL 

513:  RETURN 

514:  END 

515 : CSSSSSSSSSS  FORIGN 

516:  REAL*8  FUNCTION  FORICNCX, A.XFIX, ALP) 

517:C. . .Force  at  origin 

518:C...A  is  the  half  width 

519:  IMPLICIT  REAL*8  (A-H.O-Z) 

520:  F(T)  =  ALP*( (T-A)*DLOG(DABS(T-A))-(T+A)*DLOG(DABS(T+A)) ) 

521:  Z  =  DABS(X) 

522:  FORIGN  =  F(Z)-F(XFIX) 

523:  RETURN 

524 :  END 

525: CSSSSSSSSSS  FLNTRP 

526:  REAL*8  FUNCTION  FLNTRP(X,U,V,NDFTS,ITYP)- 

527:C. . .Function  for  linear  interpolation  on  tabular  data 
528:C. . .ITYP=1  gives  slope,  ITYP=2  gives  function  value 
529:  IMPLICIT  REAL*8  (A-H.O-Z) 

530:  DIMENSION  U(1),V(1) 

531:  IF  (X  .GT.  U(2))  GO  TO  10 

532:  K2  =  2 

533:  GO  TO  40 

534:  10  IF  (X  .LT.  U(NDPTS-l))  GO  TO  20 
535:  K2  =  NDPTS 

536:  GO  TO  40 

537:  20  K2  =  2 


538:  30  K2  =  K2  +  1 

539:  IF  (X  .GT.  U(K2))  GO  TO  30 

540:  40  K1  =  K2  -  1 

541:  IF  (DABS(U(K1 )-U(K2) )  .GT.  O.ODO)  GO  TO  50 

542:  FLNTRP  =  0.5D0  •  V(K1)  +  0.5D0  *  V(K2) 

543:  RETURN 

544:  50  IF  (ITYP.EQ.l)  FLNTRP  =  (V(K2)-V(K1) )  /  (U(K2)-U(K1 ) ) 

545:  IF  (ITYP.EQ.2) 

546:  &  FLNTRP  =  V(K1)  +  (V(K2)-V(K1))  *  (X-U(Kl))  /  (U(K2)-U(K1 ) ) 

547:  RETURN 


APPENDIX  B:  PROGRAM  II— TWO-DIMENSIONAL  STIFFNESS  MATRIX  GENERATION 
FOR  A  HALF-PLANE  OR  A  SYMMETRICALLY  LOADED  HALF-SPACE 


1.  Program  II  generalizes  the  flexibility  and  stiffness  matrix  formula¬ 
tion  used  in  Program  I.  For  problems  considering  only  vertical  loads,  stiff¬ 
ness  and  flexibility  matrices  can  be  computed  for  the  half-plane  in  plane 
strain  or  plane  stress.  An  analogous  formulation  for  a  three-dimensional  beam 
having  a  finite  width  and  resting  on  a  half-space  is  also  included.  A  more 
significant  addition  to  earlier  results  is  that  stiffness  matrices  for  the 
half-plane  subjected  to  surface  loads  involving  horizontal  components,  verti¬ 
cal  components,  and  couples  can  be  obtained.  Numerical  results  are  shown  for 
a  simple  case  involving  four  support  points,  with  the  computer  code  immedi¬ 
ately  following. 


2.  On  the  following  pages  is  a  two-dimensional  plane-strain  analysis  of 
a  sample  problem  with  vertical  surface  loads. 

Is  data  to  be  echo  printed?  (Y/N) 

N 

Is  the  program  explanation  to  be  printed?  (Y/N) 

Y 

This  program  determines  the  two-dimensional  stiffness 
matrix  for  a  series  of  loads  applied  along  the  X-axis 
on  the  surface  of  a  half-plane  or  a  half-space. 

Solutions  based  on  two-dimensional  plane  strain,  two-di¬ 
mensional  plane  stress,  or  three-dimensional  loading  can 
be  obtained.  When  a  plane  stress  or  plane -strain  condition 
is  considered,  two  types  of  loading  can  be  applied 

1)  a  series  of  vertical  loads  can  be  applied  at  a 
specified  series  of  points, 

2)  a  series  of  horizontal  loads,  vertical  loads,  and 
couples  can  be  applied. 

An  influence  function  is  used  which  computes  the 
deflection  at  any  position  due  to  a  unit  load  applied  at 
the  origin.  For  plane  strain,  a  load  of  unity  per  unit  of 
thickness  is  distributed  uniformly  between  X=-.5*A  and 
X=.5*A.  For  the  three-dimensional  loading,  a  unit  load  is 
distributed  uniformly  over  a  rectangular  zone  bounded  by 
X=-.5*A,  X=.5*A,  Z=-.5*B,  Z=.5*B  where  the  Z  direction  is 
measured  normal  to  the  XY  plane.  The  elements  of  the  stiff¬ 
ness  matrix  depend  on  E  (Young's  modulus),  NU  (Poisson's 
ratio),  the  dimension  parameters  A  and  B,  and  the  loading 
positions  X( 1 ) , . . . ,X(N) .  The  elastic  constants  enter  only 
from  a  multiplicative  factor  which  is  proportional  to  E 
for  plane  stress  or  proportional  to  E/(1.-NU**2)  for  both 
plane  strain  and  three-dimensional  loading. 

Input  values  for  Young's  modulus 
r.nd  Poisson's  ratio 
3.E5, .2 

Select  the  type  of  loading  condition 

1  =  Plane  strain 

2  =  Plane  stress 

3  =  Three  dimensional  loading 

1 

For  a  plane  strain  or  plane  stress  condition,  identify  the 
the  type  of  surface  loading 

1  =  vertical  surfaces  loads  only 

2  =  horizontal  loads,  vertical  loads,  and  couples 
l 

How  are  the  loading  positions  specified? 


r*J ■  j ’niw vr.i  FJ  V mr  V?^. vy.v.nr; r. 


1  =  Take  N  evenly  spaced  positions  between 

XMIN  and  XMAX 

2  *  Position  coordinates  X(1 ) , . . . ,X(N)  are 

input  by  the  user 

1 

Input  XMIN,  XMAX  and  the  number  of 
evenly  spaced  loading  points 
0. ,100.  ,3 

Loading  positions  are  : 

I  X(I) 

1  .0000 

2  50.0000 

3  100.0000 

For  a  plane  strain  or  plane-stress  solution 
input  the  loading  zone  width  A  (measured  in  the  X  direction) 
and  the  position  coordinate  B  at  which  zero  displacement  is  imposed 
.1,110. 

Select  compute  option: 

1  *  compute  both  flexibility  and  stiffness  matrices 

2  =  compute  flexibility . matrix  only 
1 

Select  a  print  option: 

1  =  print  flexibility  matrix  only 

2  =  print  stiffness  matrix  only 

3  =  print  both  matrices 

3 

The  flexibility  matrix  elements  on  and 
below  the  main  diagonal  are  shown  below: 

1 . 772E-05 

1 . 606E-06  1.772E-05 

1 . 942E-07  1 . 606E-06  1.772E-05 

The  stiffness  matrix  elements  on  and 
below  the  main  diagonal  are  shown  below: 

5 . 692E+04 

-5.146E+03  5 . 738E+04 

-1 . 572E+02  -5 . 146E+03  5.692E+04 


Problem  analysis  is  completed 


3.  The  solution  for  a  second  plane-strain  two-dimensional  problem  with 
horizontal  loads,  vertical  loads,  and  couples  follows. 

Is  data  to  be  echo  printed?  (Y/N) 

N 

Is  the  program  explanation  to  be  printed?  (Y/N) 

N 

Input  values  for  Young's  modulus 

and  Poisson's  ratio 

3.E5..2 

Select  the  type  of  loading  condition 

1  =  Plane  strain 

2  =  Plane  stress 

3  =  Three  dimensional  loading 

1 

For  a  plane  strain  or  plane  stress  condition,  identify  the 
the  type  of  surface  loading 

1  =  vertical  surfaces  loads  only 

2  =  horizontal  loads,  vertical  loads,  and  couples 

2 

How  are  the  loading  positions  specified? 

1  =  Take  N  evenly  spaced  positions  between 

XMIN  and  XMAX 

2  =  Position  coordinates  X(l) , . . . ,X(N)  are 

input  by  the  user 

1 

Input  XHIN,  XMAX  and  the  number  of 
evenly  spaced  loading  points 
0 . ,  100 . , 3 

Loading  positions  are  : 

I  X(I) 

1  .0000 

2  50 . 0000 

3  100 . 0000 

For  a  plane  strain  or  plane  stress  solution 
input  the  loading  zone  width  A  (measured  in  the  X  direction) 
and  the  position  coordinate  B  at  which  zero  displacement  is  imposed 
.1,110. 

Select  compute  option: 

1  =  compute  both  flexibility  and  stiffness  matrices 

2  =  compute  flexibility  matrix  only 
1 


. ■•v-v/'.y.  ■  .a 


■  «  .  «.  •  .  -  V  vo* « 


•-.A  - 


Select  a  print  option: 

1  =  print  flexibility  matrix  only 

2  =  print  stiffness  matrix  only 

3  =  print  both  matrices 
3 

The  flexibility  matrix  elements  on  and 


below  the  main  diagonal 

are  shown 

below: 

1 . 772E-05 

. 000 E +00 

1 . 772E-05 

2 . 400E-05 

.OOOE+OO 

9 . 600E-04 

1 . 606E-06 

-1 . 200E-06 

. OOOE+OO 

1 . 772E-05 

1 . 200E-06 

1 . 606E-06 

4 . 074E-08 

.OOOE+OO 

1 . 772E-05 

.OOOE+OO 

-4.074E-08 

-9 . 600E-10 

2.400E-05 

.OOOE+OO 

9 . 600E-04 

1 . 942E-07 

-1 . 200E-06 

.OOOE+OO 

1 . 606E-06 

-1 . 200E-06 

.OOOE+OO 

1 . 772E-05 

1 . 200E-06 

1 . 942E-07 

2.037E-08 

1 . 200E-06 

1 . 606E-06 

4.074E-08 

.OOOE+OO 

1 . 772E-05 

.000E+00 

-2.037E-08 

-2 . 400E-10 

.OOOE+OO 

-4 . 074E-08 

-9 . 600E-10 

2 . 400E-05 

.OOOE+OO 

9 . 600E-04 

The  stiffness  matrix  elements  on  i 

snd 

below  the  main  diagonal 

are  shown 

below: 

5.942E+04 

-1 . 237E+01 

5 . 738E+04 

-1 . 485E+03 

5 . 278E-01 

1.079E+03 

-5 . 303E+03 

3 . 739E+03 

1 . 327E+02 

5.994E+04 

-3 . 749E+03 

-4 . 937E+03 

9. 137E+01 

5.073E-14 

5. 785E+04 

1 . 327E+02 

-9 . 102E+01 

-3.319E+00 

-1.498E+03 

1.801E-15 

1.079E+03 

-4.399E+02 

3 . 324E+03 

1.084E+01 

-5. 303E+03 

3 . 749E+03 

1 . 327E+02 

5 . 942E+04 

-3 . 324E+03 

-4 . 334E+02 

8. 210E+01 

-3. 739E+03 

-4 . 937E+03 

9 . 102E+01 

1 . 237E+01 

5.738E+04 

1.084E+01 

-8 . 210E+01 

-2. 668E-01 

1 . 327E+02 

-9. 137E+01 

-3.319E+00 

-1 . 485E+03 

-5.278E-01 

1.079E+03 

Problem  analysis  is  completed 


.  The  computer  code  for  Program  II  sample  problems  is  printed  on  the 


following  pages 


CS$$$S$S$$S  MAIN 
C. . . 

C...  Program  to  determine  foundation  stiffness  matrix 

C  •  .  • 

C... Written  by:  Howard  Wilson  and  Louis  Turcotte 
C...  University  of  Alabama  at  Tuscaloosa 

C...For:  U.S.  Army  Corps  of  Engineers,  WES 

C...Date:  September  1983 

C  •  •  • 

IMPLICIT  REAL*8  (A-H.O-Z) 


LOGICAL  PRNT 

CHARACTER* 1  ANS 

COMMON  /ONE/  FLEX(60 , 60 ) , STIF( 60 , 60 ) 

COMMON  /TWO/  X ( 60 ) , STOR ( 60 ) , ISTOR ( 60 ) , FLXSTO (60,60) 

DATA  IRD/0/ , ICRT/0/ , PRNT/ .FALSE ./ , NMAX/60/ 

WRITE (ICRT, 10) 

F0RMAT(/, 

&5X, 1 •••  TWO-DIMENSIONAL  STIFFNESS  MATRIX  GENERATION  •**',/, 


&5X, ' •••  FOR  A  HALF-PLANE  OR  A  SYMMETRICALLY  •***,/, 

&5X, * •••  LOADED  HALF-SPACE  •••',/) 

WRITE (ICRT ,20) 

20  F0RMAT(/, '  Is  data  to  be  echo  printed?  (Y/N)') 

READ(IRD,30)  ANS 
30  FORMAT (Al) 

IF  (ANS  .EQ.  *Y'  .OR.  ANS  .EQ.  'y')  PRNT*. TRUE. 

IF  (PRNT)  WRITE  ( ICRT ,  AO )  ANS 
AO  FORMAT ( IX , Al ) 

WRITE ( ICRT , 50 ) 

50  FORMAT(/,'  Is  the  program  explanation  to  be  printed?  (Y/N)’) 
READ(IRD,30)  ANS 
IF  (PRNT)  WRITE (ICRT, AO)  ANS 

IF  (ANS  .NE.  'Y'  .AND.  ANS  .NE.  'y')  GO  TO  100 
WRITE (ICRT, 60) 

60  F0RMAT(/, 

&1X,'  This  program  determines  the  two-dimensional  stiffness' ,/, 
&1X, 'matrix  for  a  series  of  loads  applied  along  the  X-axis',/, 
&lX,'on  the  surface  of  a  half-plane  or  a  half-space.',/, 

&1X, ' Solutions  based  on  two-dimensional  plane  strain,  two-di-',/, 
&1X, ’ mensional  plane  stress,  or  three-dimensional  loading  can',/, 
&lX,'be  obtained.  When  a  plane  stress  or’plane  strain  condition’) 
WRITE (ICRT, 70) 

70  FORMAT ( 

&lX,'is  considered,  two  types  of  loading  can  be  applied',/, 
&1X,'1)  a  series  of  vertical  loads  can  be  applied  at  a’,/, 

&1X,'  specified  series  of  points,',/, 

&1X,'2)  a  series  of  horizontal  loads,  vertical  loads,  and',/, 
&1X,'  couples  can  be  applied.',/, 

&lX,'An  influence  function  is  used  which  computes  the') 

WRITE (ICRT, 80) 

80  FORMAT ( 

&1X, ' deflection  at  any  position  due  to  a  unit  load  applied  at',/, 


52:  &lX,’the  origin.  For  plane  strain,  a  load  of  unity  per  unit  of ' ,/, 

53:  &1X, ' thickness  is  distributed  uniformly  between  X=-.5*A  and’,/, 

54:  &1X,'X=.5*A.  For  the  three-dimensional  loading,  a  unit  load  is',/, 

55:  &1X, ' distributed  uniformly  over  a  rectangular  zone  bounded  by',/, 

56:  &1X, ' X=- . 5*A,  X=.5*A,  Z=-.5*B,  Z=.5*B  where  the  Z  direction  is',/, 

57:  &1X, ’ measured  normal  to  the  XY  plane.  The  elements  of  the  stiff-') 

58:  WRITE(ICRT,90) 

59:  90  FORMAT ( 

60:  &lX,'ness  matrix  depend  on  E  (Young''s  modulus),  NU  (Poisson' ' s' ,/ 

61:  &1X,' ratio),  the  dimension  parameters  A  and  B,  and  the  loading’,/, 

62:  &1X, 1  positions  X(l) , . . . ,X(N) .  The  elastic  constants  enter  only',/, 

63:  &lX,'from  a  multiplicative  factor  which  is  proportional  to  E' ,/, 

64:  &lX,'for  plane  stress  or  proportional  to  E/(1.-NU**2)  for  both',/, 

65:  &1X,' plane  strain  and  three-dimensional  loading.') 

66:  100  WRITE (ICRT, 110) 

67:  110  F0RMAT(/,'  Input  values  for  Young' 's  modulus',/, 

68:  &  ’  and  Poisson''s  ratio') 

69:  READ( IRD, * )  E.P0IS 

70:  IF  ( PRNT )  WRITE ( ICRT, 120 )  E,P0IS 

71:  120  FORMAT ( IX , 1PE11 . 5 , F7 . 4 ) 

72:  WRITE (ICRT, 130) 

73:  130  F0RMAT(/,*  Select  the  type  of  loading  condition’,/, 

74:  &'  1  =  Plane  strain' ,/, 

75:  &'  2  =  Plane  stress',/, 

76:  &'  3  =  Three  dimensional  loading' ) 

77:  READ (IRD,*)  ITYPE 

78:  IF  (PRNT)  WRITE(ICRT,140)  ITYPE 

79:  140  FORMAT (IX, 12) 

80:  IF  (ITYPE  .EQ.  3)  GO  TO  160 

81:  WRITE (ICRT, 150) 

82:  150  F0RMAT(/, 

83:  &'  For  a  plane  strain  or  plane  stress  condition,  identify  the',/, 

84:  &'  the  type  of  surface  loading  ' ,/, 

85:  &'  1  =  vertical  surfaces  loads  only',/, 

86:  &'  2  =  horizontal  loads,  vertical  loads,  and  couples’) 

87:  IL0AD  =  0 

88:  READ (IRD,*)  I LOAD 

89:  IF  (PRNT)  WRITE (ICRT, 140)  I LOAD 

90:  160  WRITE (ICRT, 170) 

91:  170  F0RMAT(/,'  How  are  the  loading  positions  specified?',/, 

92:  &'  1  =  Take  N  evenly  spaced  positions  between’ ,/, 

93:  &'  XMIN  and  XMAX’ ,/, 

94:  &'  2  =  Position  coordinates  X( 1 ) , . . . ,X(N)  are',/, 

95:  &'  input  by  the  user') 

96:  READ (IRD,*)  I0PT 

97:  IF  (PRNT)  WRITE (ICRT, 140)  10 PT 

98:  IF  ( IOPT  .EQ.  2)  GO  TO  240 

99:  WRITE ( ICRT, 180) 

100:  180  FORMAT ( / , '  Input  XMIN,  XMAX  and  the  number  of’,/, 

101:  &'  evenly  spaced  loading  points') 

102:  READdRD,*)  XMIN, XMAX, N 

103:  IF  (PRNT)  WRITE ( ICRT, 190)  XMIN, XMAX, N 

104:  190  FORMAT (IX, F10. 4, ' , ' ,F10.4, ' , ’ ,13) 

105:  D  =  (XMAX -XMIN)  /  (N-l) 
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DO  200  1=1, N 

200  X(I)  =  XMIN  +  D  *  (1-1) 

WRITE (ICRT, 210) 

210  F0RMAT(/, '  Loading  positions  are  I* ,5X, 'X(I) ' ) 

DO  220  1=1, N 

220  WRITE (ICRT, 230)  I,X(I) 

230  FORMAT ( IX , 13 , F10 , 4 ) 

GO  TO  270 

240  WRITE (ICRT, 250) 

250  F0RMAT(/,'  Input  N  and  X(l) . X(N>') 

READdRD,*)  N,(X(K)  ,K=1,N) 

IF  (PRNT)  WRITE( ICRT, 260)  N, (X(K) ,K=1,N) 

260  F0RMAT(1X,I3,/, 

&  7 ( F10, 4 , ' , ' , F10. 4 , ' , ' , F10.4 , * , ‘ , F10. 4 ,  ’  , ' ,F10. 4 , ' , ' ,F10. 4 ,/) ) 

270  IF(ITYPE.EQ.3)  GO  TO  300 
WRITE (ICRT, 280) 

280  F0RMAT(/,'  For  a  plane  3train  or  plane  stress  solution',/, 

&'  input  the  loading  zone  width  A  (measured  in  the  X  direction)', 
&/,'  and  the  position  coordinate  B  at  which  zero  displacement', 

&’  is  imposed') 

READdRD,*)  A,B 

IF  (PRNT)  WRITE (ICRT, 290)  A,B 

290  FORMAT ( IX , F10 . 4 , ' , ' , F10 . 4 ) 

GO  TO  315 

300  WRITE (ICRT, 310) 

310  F0RMAT(/,'  For  the  dimensions  of  loading  input  the',/, 

&'  X-width  and  the  Z-width  of  the  loading  zone') 

READdRD,*)  A,B 

IF  (PRNT)  WRITE (ICRT, 2 90)  A,B 

315  WRITE (ICRT, 320) 

320  F0RMAT(/,'  Select  compute  option:',/, 

&  '  1  =  compute  both  flexibility  and  stiffness  matrices',/, 

&  '  2  =  compute  flexibility  matrix  only' ) 

READdRD,*)  NOSTIF 

IF  (PRNT)  WRITE (ICRT, 140)  NOSTIF 

NOSTIF  =  NOSTIF  -  1 

330  IF  (ILOAD  .EQ.  1  .OR.  ITYPE  .EQ.  3  )  THEN 

CALL  STFVLO ( X , N , NMAX , E , POIS , A , B , ITYPE , STOR , ISTOR , IERROR , NOSTIF , 

&  FLXSTO , FLEX , STI F ) 

ELSE 

CALL  STFHVM ( X , N , NMAX , E , POIS , A , B , ITYPE , STOR , ISTOR , IERROR , NOSTIF , 

&  FLXSTO, FLEX, STI F) 

ENDIF 

IF  (IERROR  .NE.  1)  GO  TO  350 
WRITE (ICRT, 340) 

340  F0RMAT(/, '  A  SINGULAR  FLEXIBILITY  MATRIX  WAS’,/, 

&’  OBTAINED.  EXECUTION  IS  TERMINATED',/) 

STOP  ’  ’ 

350  WRITE(ICRT,360) 

360  FORMAT(/,’  Select  a  print  option:’, 

&  /,’  1  =  print  flexibility  matrix  only', 

&  /,'  2  =  print  stiffness  matrix  only', 

&  /,'  3  =  print  both  matrices') 

READdRD,*)  NOPT 


B8 


160:  IF  (PRNT)  WRITE( ICRT, 140)  NOPT 

161:  NN  =  N 

162:  IF  (ILOAD  .EQ.  2)  NN  =  3  *  N 

163:  IF  (NOPT  .EQ.  2)  GO  TO  390 

164 :  WRITE (ICRT, 370) 

165:  370  FORMAT^/,'  The  flexibility  matrix  elements  on  and’,/, 

166:  &’  below  the  main  diagonal  are  shown  below:') 

167:  DO  380  1*1, NN 

168:  380  WRITE (ICRT, 420)  (FLEX(I,K) ,K=1,I) 

169:  390  IF  (NOPT  .EQ.  1)  GO  TO  430 
170:  WRITE ( ICRT, 400) 

171:  400  FORMAT(/,'  The  stiffness  matrix  elements  on  and  ' ,/, 

172:  &'  below  the  main  diagonal  are  shown  below:') 

173:  DO  410  1=1, NN 

174:  410  WRITE (ICRT, 420)  (STIF(I,K) ,K=1,I) 

175:  420  F0RMAT(1X,6(1PE11.3) ,/,12X,5(lPE11.3) ,/,12X,5(lPE11.3) ,/, 

176:  &  12X,5(1PE11.3) ,/,12X,5(lPE11.3) ,/ , 12X,5 ( 1PE11 . 3 ) ,/, 

177:  &  12X, 5(1PE11.3) ,/ , 12X, 5 (1PE11 . 3) ,/,12X, 5(1PE11.3) ) 

178:  430  WRITE (ICRT, 440) 

179:  440  F0RMAT(/,'  Problem  analysis  is  completed',/) 

180:  STOP  '  ’  . 

181:  END 

182:C$S$$$$$$$S  STFVLO 

1 83 :  SUBROUTINE  STFVLO (X,N, NMAX , E , POIS , A , B , ITYPE , STOR , ISTOR , NOSTIF , 

184:  &  I ERROR , FLXSTO , F  LEX , ST I F ) 

185:C. . .Stiffness  ma5rixrand  flexibility  formulation  for  half-plane 
186:C...or  half-space  subjected  to  vertical  loads  only. 

187 :C. . . 

188 :C. . .E  =  Young's  modulus 

189 :C... POIS  =  Poisson's  ratio. 

190:C...A,  B  =  the  loading  dimension  parameters 

191:C...X(1),...,X(N)  are  the  load  application  positions 

192 : C. .. ITYPE  =  1,2,  or  3  for  plane  strain,  plane  stress,  or 

193:C...  three-dimensional  loading 

194 :C. . .NOSTIF  =  1  to  generate  only  flexibility  matrix 

195:C...  0  to  generate  both  flexibility  and  stiffness  matrix 

196:C...FLEX  =  the  flexibility  matrix 

197:C...STIF  =  the  stiffness  matrix  obtained  by  inverting  the  flexibility 

198 :C...  matrix  when  NOSTIF  *  0 

199 :C. . .IERROR  =  1  for  a  normal  return 

200 :C...  2  for  flexibility  matrix  that  is  singular 

201 :C... STOR  =  working  storage  of  length  N  or  more 

202 : C ... ISTOR  =  working  storage  of  length  N  or  more 

203:C...STIF  =  the  N  by  N  roundation  stiffness  matrix  which  is  the 

204 :C...  desired  subroutine  output 

205:  IMPLICIT  REAL*8  (A-H.O-Z) 

206:  DIMENSION  X(N) ,FLEX(NMAX,N) ,STIF(NMAX,N) ,ST0R(1 ) ,IST0R(1 ) , 

207:  &  FLXSTO (NMAX,N) 

208:  DATA  ICRT/  0  / 

209:  DO  20  1*1, N 

210:  XI  *  X(I) 

211:  DO  10  J*1,I 

212:  XIJ  *  XI  -  X(J) 

213:  F  *  FUNIT(XIJ, A, B,E, POIS, ITYPE) 


214:  FLEX(I.J)  *  F 

215:  FLEX(J.I)  *  F 

216:  FLXSTO(I.J)  =  F 

217:  10  FLXSTO(J.I)  *  F 

218:  20  CONTINUE 

219 :C. . .Invert  the  flexibility  matrix  to  get  the  stiffness  matrix 
220:  IERR0R-0 

221:  CALL  INVERT ( FLXSTO , N , NMAX , ISTOR , STOR , IFLAG , STIF ) 

222:C...Make  certain  stiffness  matrix  is  symmetric 
223:  DO  40  1*1, N 

224:  DO  40  J*1,I 

225:  SIJ  *  0.5D0  •  (  STIF(I.J)  +  STIF(J.I)  ) 

226:  STIF(I.J)  *  SIJ 

227:  40  STIF(J.I)  *  SIJ 

228:C. . .Check  whether  return  from  INVERT  was  normal 

229:  IF  (IFLAG  .EQ.  2)  I ERROR  -  1 

230:  RETURN 

231 :  END 

232 : CSSSSSSSSSS  FUNIT 

233:  REAL*8  FUNCTION  FUNIT(X,A,B,E,POIS, ITYPE) 

234 :C. . .This  function  returns  the  deflection  at  position  X  due 
235:C...to  a  unit  distributed  load  at  the  origin.  Parameter 
236 : C . . . ITYPE  equals  1,  2,  or  3  depending  on  whether  a  plane 
237:C. . .strain,  a  plane  stress,  or  a  three-dimensional  solution 
238:C...is  generated.  E  and  POIS  denote  Young's  modulus  and 
239:C. . .Poisson' s  ratio,  respectively.  In  instances  of  plane 
240 : C. . .strain  or  plane  stress,  the  unit  load  is  distributed 
241;C. . .over  an  area  of  width  A  and  unit  depth.  The  displace- 
242:C...ment  is  adjusted  to  equal  zero  at  X*B  in  the  plane 
243:C. . .case.  In  the  instance  of  three-dimensional  loading, 

244 :C... the  unit  load  is  distributed  over  a  rectangular  area 
245:C...of  width  A  and  depth  B.  The  displacement  at  X*infinity 
246:C. . .vanishes  for  the  three-dimensional  case. 

247:  IMPLICIT  REAL*8  (A-H.O-Z) 

248:  DATA  PI/  3. 14159265358979DO  / 

249:  IF  (ITYPE  .EQ.  3)  GO  TO  10 

250 :C. . .Plane  stress 

251:  C  *  2. DO  /  (  PI  *  E  *  A  ) 

252 :C. . .Plane  strain 

253:  IF  (ITYPE  .EQ.  1)  C  -  C  •  (  1.D0  -  P0IS*P0IS  ) 

254:  FUNIT  *  -C  *  (F2D(DABS(X) ,.5D0*A)-F2D(B, .5D0*A)) 

255 :  RETURN 

256 : C. . .Three  dimensional  case 

257:  10  T2  *=  (  A  +  2. DO  *  DABS(X)  )  /  B 

258:  T1  ■  (  -A  +  2. DO  *  DABS(X)  )  /  B 

259:  20  FUNIT  -  (F3D(T2)-F3D(T1))  •  (1.D0-P0IS*P0IS)  /  (PI*E*A) 

260:  RETURN 

261:  END 

262:C$$$$$$$$$$  F2D 

263:  REAL*8  FUNCTION  F2D(U,V) 

264:  IMPLICIT  REAL*8  (A-H.O-Z) 

265:  UA  *  DABS(U) 

266:  F2D  -  (UA+V)  •  DL0G(UA+V) 

267:  IF  (UA  .NE.  V)  F2D  *  F2D  -  (UA-V)  •  DLOC(DABS(UA-V)) 


268:  RETURN 

269:  END 

270:C$$$$$$S$$$  F3D 

271:  REAL*8  FUNCTION  F3D(T) 

272:  IMPLICIT  REAL*8  (A-H.O-Z) 

273:  F3D  =  DLOG(T+DSQRT(l .DO+T*T) ) 

274:  IF  (T  .EQ.  O.DO)  RETURN 

275:  S  =  1 . DO  /  DABS(T) 

276:  F3D  =  F3D  +  T  *  DLOG(S+DSQRT(1.DO+S*S) ) 

277:  RETURN 

278:  END 

279 : CSSSSSSSSSS  FACTR 

280:  SUBROUTINE  FACTR  (A, N, NMAX, IPIVOT , SCAL, IFLAG) 

281 : C. .. Perform  triangular  factorization  of  matrix  A 
282: C. . .using  scaled  row  pivoting. 

283:C. .. IFLAG  =  1  means  normal  return 
284:C...  =  2  means  matrix  is  singular 

285:  IMPLICIT  REAL*8  (A-H.O-Z) 

286:  DIMENSION  A(NMAX,N) ,  IPIVOT(N),  SCAL(N) 

287:  IFLAG  =  1 

288:C. . .Initialize  IPIVOT  and  SCAL 

289:  DO  20  1=1. N 

290:  IPIVOT(I)  =  I 

291:  ROWMAX  =  O.DO 

292:  DO  10  J=1,N 

293:  10  ROWMAX  =  DMAX1 (ROWMAX, DABS (A( I ,J) ) ) 

294:  IF  (ROWMAX  .EQ.  O.DO)  GO  TO  50 

295:  20  SCAL(I)  =  ROWMAX 

296:C. .. Perform  Gauss  reduction 

297:  NM1  =  N  -  1 

298:  IF  ( NM1  .EQ.  0)  RETURN 

299:  DO  40  K=1,NM1 

300:  J  =  K 

301:  KP1  =  K  +  1 

302:  IP  =  IPIVOT(K) 

303:  COLMAX  =  DABS(A(IP,K) )  /  SCAL(IP) 

304:  DO  30  I=KP1,N 

305:  IP  =  IPIVOT ( I ) 

306:  AWIKOV  =  DABS( A( IP, K) )  /  SCAL(IP) 

307:  IF  (AWIKOV  .LE.  COLMAX)  GO  TO  30 

308:  COLMAX  =  AWIKOV 

309:  J  =  I 

310:  30  CONTINUE 

3x1:  IF  (COLMAX  .EQ.  O.DO)  GO  TO  50 

312 ;  C 

313:  IPK  =  IPIVOT ( J ) 

314:  IPIVOT ( J )  =  IPIVOT (K) 

315:  IPIVOT(K)  =  IPK 

316:  DO  40  '=KP1,N 

317:  IP  =  IPIVOT ( I ) 

318:  A( IP,K)  =  A( IP, K)  /  A(IPK,K) 

319:  RATIO  =  -A(IP,K) 

320:  DO  40  J=KP1,N 

321:  40  A( IP , J )  =  RATIO  •  A(IPK,J)  +  A(IP,J) 
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322:  IF  (A(IP,N)  .EQ.  O.DO)  GO  TO  50 

323:  RETURN 

324 :  50  IFLAG  =  2 

325:  RETURN 

326 :  END 

327 :C$$$$SS$$$$  BAKSUB 

328:  SUBROUTINE  BAKSUB  (A,N,NMAX,B, IPIVOT, X) 

329:C. . .Solve  simultaneous  equations  AX=B  where  matrix  A  has 
33Q:C...been  subjected  to  factorization  by  subroutine  FACTR 
331:  IMPLICIT  REAL*8  (A-H.O-Z) 

332:  DIMENSION  A(NMAX,N) ,  B(N) ,  IPIVOT(N),  X(N) 

333:  IF  (N  .GT.  1)  GO  TO  10 

334:  X(l)  a  B(l)  /  A(l,l) 

335 :  RETURN 

336:  10  IP  =  IPIVOT(l) 

337:  X(l)  =  B(IP) 

338:  DO  30  K*2,N 

339:  IP  a  IPIVOT(K) 

340:  KM1  *  K  -  1 

341:  SUM  a  O.DO 

342:  DO  20  J-l.KMl 

343:  20  SUM  =  A(IP,J)  *  X(J)  +  SUM 

344:  30  X(K)  =  B(IP)  -  SUM 

345:  X(N)  a  X(N)  /  A(IP,N) 

346:  K  a  N 

347;  DO  50  NP1MK=2,N 

348:  KP1  =«  K 

349:  K  *  K  -  1 

350:  IP  «  IPIVOT(K) 

351:  SUM  =  O.DO 

352:  DO  40  J=KP1,N 

353:  40  SUM  a  A(IP,J)  *  X(J)  +  SUM 

354:  50  X(K)  a  (X(K)-SUM)  /  A(IP,K) 

355 :  RETURN 

356:  END 

357 :C$$$$$$$$$$  INVERT 

358:  SUBROUTINE  INVERT (  A,  N,  NMAX,  IPIVOT,  SCAL,  IFUG,  AINV  ) 

359 :C. .. IFLAG  a  2  indicates  singular  matrix 
360:  IMPLICIT  REAL*8  (A-H,0-Z) 

361:  DIMENSION  A(NMAX,1),  AINV (NMAX, 1) ,  IPIVOT(l),  SCAL(l) 

362:  CALL  FACTR (  A,  N,  NMAX,  IPIVOT,  SCAL,  IFUG  ) 

363:  IF (  IFUG  .EQ.  2  )  RETURN 

364:  DO  10  J  a  1,  N 

365:  10  SCAL(J)  =  O.ODO 

366:  DO  20  I  a  1,  N 

367:  SCAL(I)  a  l.ODO 

368:  CALL  BAKSUB (  A,  N,  NMAX,  SCAL,  IPIVOT,  AINV(1,I)  ) 

369:  SCAL ( I )  =  O.ODO 

370 :  20  CONTINUE 

371 :  RETURN 

372:  END 

373:C$S$$S$$$$$  F 

374:  REAL*8  FUNCTION  F(  X,  A  ) 

375 :C... Integral  from  0  to  X  of  LOG(ABS(X-A) )-LOG(ABS(X+A) ) 
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376:  IMPLICIT  REAL*8  (A-H.O-Z) 

377:  FF  =  O.DO 

378:  XMA  =  DABS(X)  -  A 

379:  XPA  =  DABS(X)  +  A 

380:  IF  (XMA  .E Q.  O.DO)  GO  TO  10 

381:  FF  =  XMA  *  DLOG (DABS (XMA ) ) 

382:  10  FF  =  FF  -  XPA  *  DLOG(XPA) 

383:  F  =  FF 

384:  RETURN 

385 :  END 

386 :C$$$$$$$$$$  STEP 

387:  REAL*8  FUNCTION  STEP(  X,  A  ) 

388:C...Unit  step  function. 

389:  IMPLICIT  REAL*8  (A-H.O-Z) 

390:  STEP  =  O.DO 

391:  IF  (X  .GE.  A)  STEP  =  l.DO 

392:  RETURN 

393 :  END 

394 :C$$$$$$$S$$  RAMP 

395:  REAL*8  FUNCTION  RAMP(  X  ,  A  ) 

396 :C. .. Linearly  varying  ramp  load. 

397:  IMPLICIT  REAL*8  (A-H.O-Z) 

398:  RAMP  =  O.DO 

399:  IF  (X  .GE.  A)  RAMP  =  X  -  A 

400 :  RETURN 

401 :  END 

402 : C$S$$$S$$$S  R 

403:  REAL*8  FUNCTION  R(  X,  A  ) 

404 :C. . .Derivative  of  F(X,A) 

405:  IMPLICIT  REAL*8  (A-H.O-Z) 

406:  DATA  ICKT  /  0  / 

407:  XMA  =  DABS(X-A) 

408:  IF  (XMA  .EQ.  O.DO)  GO  TO  10 

409:  R  =  DLOG (XMA) 

410:  XPA  =  DABS(X+A) 

411:  IF  (XPA  .EQ.  O.DO)  GO  TO  10 

412:  R  =  R  -  DLOG(XPA) 

413:  RETURN 

414:  10  WRITE (  ICRT,  20  ) 

415:  20  FORMAT (  /,  '  ARGUMENT  ERROR  IN  FUNCTION  RAMP',  /  ) 

416:  STOP  '  ' 

417:  END 

418:C$$$$$$$$$$  S 

419:  REAL *8  FUNCTION  S(  X,  A  ) 

420:C. .. Second  derivative  of  F(X,A) 

421:  IMPLICIT  REAL*8  (A-H.O-Z) 

422:  DATA  ICRT  /  0  / 

423:  IF  (DABS(X)  .NE.  DABS(A))  GO  TO  20 

424:  WRITE (  ICRT,  10  ) 

425;  10  FORMAT (  /,  *  ARGUMENT  ERROR  IN  FUNCTION  S’,  /  ) 

426:  STOP  ’  ' 

427:  20  S=2.D0*A/(X*X-A*A) 

428:  RETURN 

429 


430:0$$$$$$$$$$  G 

431:  REAL*8  FUNCTION  G(  X,  A  ) 

432 :C... Integral  of  H(X,A) 

433:  IMPLICIT  REAL*8  (A-H.O-Z) 

434:  G  =  -A 

435:  IF  (X  .GT.  -A)  G  =  X 

436:  IF  (X  .GT.  A)  G  -  A 

437:  RETURN 

438:  END 

439 : C $$$$$$$$$$  H 

440:  REAL*8  FUNCTION  H(  X,  A  ) 

441 :C. . .Function  equals  1  between  -A  and  +A  and  equals  zero  otherwise. 

442:  IMPLICIT  REAL*8  (A-H.O-Z) 

443:  H  =  O.DO 

444:  IF  (DABS(X)  .LE.  A)  H  =  1.D0 

445:  RETURN 

446:  END 

447 :C$$$$$$S$$$  STFHVM 

448:  SUBROUTINE  STFHVM(X, N.N3MAX, E, POIS, WIDTH, DFIX, ITYP, STOR, ISTOR, 

449:  &  IERROR,NOSTIF,FLXSTO,FLEX,STIF) 

450:C...This  subroutine  generates  flexibility  and  stiffness  matrices  for 
451:C...a  half  plane  subjected  to  horizontal  loads,  vertical  loads,  and 
452:C. . .couples  applied  at  a  selected  set  of  coordinates  on  the  surface. 
453:C. . . 

454:C...This  routine  calls  the  following  routines:  F,  STEP,  RAMP,  R,  S, 
455:C...C,  H,  FACTR,  BAKSUB,  INVERT 
456:  IMPLICIT  REAL*8  (A-H.O-Z) 

457:  DIMENSION  X(N),  FLEX ( N 3M AX , N3MAX ) ,  ST0R(1),  ISTOR(l), 

458:  &  FLXSTO ( N 3MAX , N 3MAX > ,  STIF(N3MAX,N3MAX) 

459:  DATA  PI  /  3. 14159265358979D0  /,  ICRT/  0  / 

460:  AH  *  0.5D0  *  WIDTH 

461:  FD  =  F(DFIX, AH) 

462:  IF  (ITYP  .EQ.  2)  GO  TO  10 

463 :C. .. Plane  strain 

464:  CO  =  (  1 . DO  -  POIS  **  2  )  /  (  PI  *  AH  •  E  ) 

465:  EO  =  (  1.D0  +  POIS  )  *  (  1.D0  -  2. DO  •  POIS  )  /  (  2. DO  *  AH  *  E  ) 

466:  GO  TO  20 

467:C. . .Plane  stress 

468:  10  CO  a  l.DO  /  (  E  *  PI  *  AH  ) 

469:  EO  =  (  l.DO  -  POIS  )  /  (  2. DO  *  E  •  AH  ) 

470:  20  CONTINUE 

471:  DO  40  I  a  1,  N 

472:  DO  30  J  =  1,  N 

473:  XIJ  =  X(I)  -  X(J) 

474:  FIJ  »  CO  *  (  F(XIJ,AH)  -  FD  ) 

475:  GIJ  =  EO  *  G(  XIJ,  AH  ) 

476:  HIJ  -  EO  *  H(  XIJ,  AH  ) 

477:  RIJ  a  CO  •  R(  XIJ,  AH  ) 

478:  SIJ  *  EO  *  S(  XIJ,  AH  ) 

479:  IR  *  3  •  I  -  2 

480:  JC  «  3  *  J  -  2 

481:  FLEX(IR.JC)  *  FIJ 

482:  FLEX(IR+1 , JC)  =  GIJ 

483:  FLEX(IR+2, JC)  =  HIJ 
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484:  FLEX(IR, JC+1 )  =  -GIJ 

485:  FLEX (IR+1, JC+1)  =  FIJ 

486:  FLEX (IR+2, JC+1)  =  RIJ 

487:  FLEX(IR, JC+2)  =  HIJ 

488:  FLEX (IR+1, JC+2)  =  -RIJ 

489:  FLEX (IR+2, JC+2)  =  -SIJ 

490:  30  CONTINUE 

491:  40  CONTINUE 

492:  N3  =  3  *  N 

493:  DO  50  1=1, N3 

494:  DO  50  J=1,I 

495:  FSTO  =  FLEX(I,J) 

496:  FLXSTOd,  J)  =  FSTO 

497:  50  FLXSTOd,  I)  =  FSTO 

498 :C. .. Invert  the  flexibility  matrix  if  stiffness  is  required 
499:  IF  (NOSTIF  .EQ.  1)  RETURN 

500:  I ERROR  =  0 

501:  CALL  INVERT ( FLXSTO , N3 , N3MAX , ISTOR , STOR , IFLAG, STIF ) 

502:C...Make  sure  stiffness  matrix  is  symmetric 
503:  DO  60  1=1, N3 

504:  DO  60  J=1,I 

505:  SIJ  =  0.5D0  *  (  STIF(I,J)  +  STIF(J.I)  ) 

506:  STIF(I,J)  =  SIJ 

507:  60  STIF(J.I)  =  SIJ 

508 : C. . .Check  whether  return  form  INVERT  was  normal 
509:  IF  (IFLAG  .EQ.  2)  I ERROR  =  1 

510:  RETURN 

511:  END 


